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Abstract 


Parameter  estimation  is  a  key  component  of  many  flight 
test  programs,  particularly  during  the  initial  development 
of  new  designs.  Accurate  estimates  of  key  aerodynamic 
derivatives  are  required  to  update  wind  tunnel  models  and 
simulations.  These  estimates  facilitate  envelope  expansion 
flights,  flight  control  system  design,  and  increase  fidelity 
of  ground  simulators.  With  the  advent  of  more  highly 
augmented  aircraft,  more  sophisticated  techniques  are 
required  (over  and  above  conventional  analysis)  to  provide 
accurate  estimates  of  aerodynamic  derivatives.  While  there 
are  computer  programs  which  have  been  used  for  years  at  the 
Air  Force  Flight  Test  Center  (AFFTC) ,  such  as  Modified 
Maximum  Likelihood  Estimation  (MMLE3),  they  are  mainframe 
computer  based,  and  require  considerable  time  to  set  up  for 
a  particular  test  program. 

This  thesis  develops  and  tests  a  personal  computer 
based  approach  to  the  aircraft  parameter  identification 
problem.  The  approach  was  tested  analytically  at  the  Air 
Force  Institute  of  Technology,  Wright-Patterson  AFB,  Ohio. 
These  tests  were  conducted  using  simplified  short  period 
approximation  and  second-order  transfer  function  models  of 
longitudinal  aircraft  response.  Initial  evaluations  were 
accomplished  using  simulated,  noise-free  data. 

xvi 


The  longitudinal  aircraft  model  was  modified  prior  to 
flight  test,  for  both  convenience  and  necessity.  This  more 
complex  model  was  evaluated  during  a  flight  test  program  at 
the  USAFTPS,  Edwards  AFB,  California.  A  T-38A  aircraft  was 
used  as  a  tool  to  demonstrate  and  validate  the  process. 
Stability  derivative  estimates  from  the  personal  computer 
based  approach  were  compared  with  those  from  the  mainframe 
MMLE3  program.  Second-order  system  parameters  calculated 
using  these  derivatives  were  also  compared  to  those  obtained 
from  frequency  response  analysis  routines  also  available  on 
the  Cyber  computer  at  the  AFFTC. 
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DEVELOPMENT  OF  A  PERSONAL  COMPUTER  BASED  APPROACH 


TO  AIRCRAFT  PARAMETER  IDENTIFICATION 


I .  Introduction 


Background 

Parameter  estimation  is  a  key  component  of  many  flight 
test  programs,  particularly  during  the  initial  development 
of  new  designs.  Accurate  estimates  of  key  aerodynamic 
derivatives  are  required  to  update  wind  tunnel  models  and 
simulations.  These  estimates  facilitate  envelope  expansion 
flights,  flight  control  system  design,  and  increase  fidelity 
of  ground  simulators.  With  the  advent  of  more  highly 
augmented  aircraft,  more  sophisticated  techniques  are 
required  (over  and  above  conventional  analysis)  to  provide 
accurate  estimates  of  aerodynamic  derivatives. 

Conventional  dynamic  techniques  consist  of  maneuvers 
such  as  stick  and  rudder  pedal  pulses  or  doublets.  The  time 
history  responses  of  the  aircraft  are  analyzed  using  log 
decrement  or  time  ratio  methods  to  determine  damping  ratio 
and  natural  frequency  (1:8.47-8.53).  In  general,  these 
conventional  techniques  have  three  primary  problem  areas: 
accuracy,  ease,  and  speed  of  testing  and  analysis. 
Conventional  testing  requires  specific  maneuvers  for  each 
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parameter  estimate.  In  some  cases,  highly  augmented 
aircraft  have  control  systems  which  prevent  analysis  of 
basic  aircraft  motion.  Even  in  those  cases  where  basic 
aircraft  motion  can  be  analyzed,  the  process  is  time 
consuming  and  often  inaccurate.  Finally,  from  the 
standpoint  of  flight  safety,  highly  augmented  control 
systems  could  mask  serious  handling  qualities  problems, 
making  it  difficult  (and  perhaps  dangerous)  to  extrapolate 
handling  qualities  trends  (11:11). 

Derivative  analysis  techniques  offer  several  advantages 
over  conventional  techniques.  First  of  all,  derivatives 
provide  better  insight  into  handling  qualities  discrepancies 
between  predicted  and  actual  flight  test  data.  Because  of 
the  dependence  of  many  conventional  flight  test 
relationships  on  multiple  stability  derivatives,  it  is 
difficult  to  determine  which  derivative  is  in  error  if  a 
discrepancy  occurs.  Derivative  analysis  also  allows  for 
extrapolation  and  standardization  of  data.  Extrapolation 
allows  for  a  safer  approach  to  envelope  expansion  and  high 
angle  of  attack  testing.  Standardization  makes  it  possible 
to  make  direct  comparisons  of  data  from  different  tests  or 
different  aircraft.  Techniques  for  derivative  analysis  can 
be  easily  managed  in  computers,  and  the  results  obtained  can 
be  used  to  update  simulators  with  flight  test  data. 

Finally,  there  is  an  overall  savings  in  flight  test  time. 
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This  is  becoming  especially  important  as  costs  associated 
with  flight  test  increase  (11:19-21). 

Objective 

The  lack  of  a  flexible  implementation  of  parameter 
estimation  has  been  a  stumbling  block  in  the  past.  While 
there  are  computer  programs  which  have  been  used  for  years 
at  the  Air  Force  Flight  Test  Center  (AFFTC) ,  such  as 
Modified  Maximum  Likelihood  Estimation  (MMLE3) ,  they  are 
mainframe  computer  based,  and  require  considerable  time  to 
set  up  for  a  particular  test  program.  Recently,  a 
commercially  available  personal  computer  software  package 
has  become  available  which  implements  the  MMLE3  problem 
formulation  within  the  PC  Matlab  program.  PC  Matlab  is  a 
matrix  operation  program  which  can  be  used  to  solve 
practical  engineering  and  mathematical  problems.  PC  Matlab 
is  designed  for  specific  applications  such  as  control  theory 
and  signal  processing.  (8:3).  Two  additional  "toolboxes" 
are  used  for  this  project:  the  Control  System  Toolbox  (6) 
and  the  State  Space  Identification  Tool  (10).  The  State 
Space  Identification  Tool  implements  the  MMLE3  formulation, 
but  it  is  written  in  terms  of  general  system  identification, 
and  is  not  designed  for  specific  application  to  aircraft 
equations  of  motion.  The  objective  of  this  thesis  is  to 
modify  and  test  the  existing  personal  computer  software 
package  to  accomplish  aircraft  parameter  estimation.  The 
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end  goal  is  to  verify  and  validate  the  results  of  the 
personal  computer  based  MMLE  by  comparing  flight  test 
results  with  those  from  the  mainframe  MMLE3  program,  as  well 
as  cross  checking  the  results  with  conventional  analysis 
techniques  and  frequency  response  analysis  methods. 

Section  II  will  discuss  the  theory  behind  parameter 
estimation,  including  maximum  likelihood  estimation  and  some 
of  the  particular  capabilities  of  the  existing  Fortran 
program,  MMLE3.  Section  III  presents  the  work  done  using 
simulated  data  to  check  out  the  PC  Matlab  MMLE  software. 
Section  IV  summarizes  the  final  preparations  which  were 
necessary  prior  to  proceeding  with  flight  test.  Section  V 
describes  the  details  of  the  flight  test  program,  which  was 
the  ultimate  test  of  the  personal  computer  based  approach. 
Finally,  section  VI  summarizes  the  conclusions  drawn  and 
recommendations  made  as  a  result  of  this  effort. 


4 


II.  Parameter  Estimation  Theory 


The  general  problem  of  aircraft  parameter  estimation 
can  be  described  in  fairly  simple  terms.  It  is  first 
assumed  that  a  set  of  dynamic  equations  exists  which 
accurately  model  the  system  in  question.  This  model 
contains  various  unknowns,  which  are  the  parameters  of 
interest.  The  system  is  then  subjected  to  a  predefined 
input,  designed  to  produce  a  response  which  will  allow 
estimation  of  these  unknowns.  Once  the  actual  input  and 
responses  are  measured,  values  for  the  unknown  parameters 
are  derived  from  the  constraint  that  the  model  responses 
should  match  the  actual  responses.  This  whole  process 
sounds  fairly  simple,  but  several  problems  arise  when  trying 
to  apply  the  theory  to  a  real  system. 

This  first  problem  is  measurement  noise.  It  is 
basically  impossible  to  measure  a  real  system  perfectly. 
Measurement  noise  makes  it  impossible  to  identify  actual 
values  of  the  parameters  of  interest.  This  causes  the 
problem  to  be  considered  parameter  estimation,  rather  than 
parameter  identification. 

The  second  problem  is  state  noise.  State  noise  is  any 
outside  (unmeasured)  source  that  causes  unexpected  system 
response.  For  the  aircraft  problem,  atmospheric  turbulence 
is  the  most  common  form  of  state  noise.  With  both 
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measurement  and  state  noise,  the  problem  becomes  more 
difficult  than  with  either  of  the  two  alone. 

Finally,  modeling  error  can  corrupt  the  reliability  of 
the  parameter  estimates.  In  basic  terms,  it  is  assumed 
throughout  the  parameter  estimation  process  that  there  is  a 
dynamic  model  which  accurately  represents  the  system 
response.  In  other  words,  for  some  set  of  "true”  values  for 
each  unknown  parameter,  it  is  assumed  that  the  model  and 
actual  responses  match.  In  most  cases,  simplified  dynamic 
models  do  not  exactly  describe  the  real  system. 

For  purposes  of  aircraft  parameter  identification,  the 
MMLE3  program  has  been  used  for  many  years  as  the  standard 
at  the  AFFTC.  The  basics  of  the  MMLE3  program  lie  in  the 
theory  of  Maximum  Likelihood  Estimation,  which  is  reviewed 
in  the  following  subsection. 

Maximum  Likelihood  Estimation 

The  Modified  Maximum  Likelihood  Estimation  {MMLE3) 
program  is  a  Fortran  computer  program  which  is  used  as  the 
standard  for  stability  derivative  extraction  at  the  AFFTC. 
This  subsection  describes  the  method  used  by  MMLE3  for 
maximum  likelihood  estimation,  as  developed  by  Maine  and 
Iliff  (7) .  In  general,  the  parameter  estimates  are  obtained 
by  choosing  a  parameter  vector  f  which  maximizes  the 
likelihood  functional  p(z|f)»  where  z  is  the  measured  system 
response.  The  critical  part  of  the  development  is  the 
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actual  definition  of  p(z|f),  which  leads  directly  to  the 
cost  function  for  the  minimization  process.  With  both  state 
and  measurement  noise,  the  continuous  system  model  is 

x(t)  =  Ax(t)  +  Su(t)  +  Fn(t) 


where 


X  =  state  vector 
u  =  input  vector 
n  =  state  noise  vector 
A  =  continuous  system  dynamics  matrix 
B  =  control  distribution  matrix 

F  =  square  root  of  the  state  noise  spectral  density 


The  discrete  measurement  equation  is 

=  Cx(t^)  +  Du(tj)  +Gti(tj) 


(2) 


where 

z  =  measured  observation  vector 
rj  =  measurement  noise  vector 
C  =  state  measurement  matrix 
D  =  discretized  transmission  matrix 
G  =  square  root  of  the  measurement  noise  covariance 
matrix 
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It  is  assumed  that  the  state  noise  vector  is  zero  mean, 
Gaussian,  and  white,  with  unit  spectral  density.  Also,  the 
measurement  noise  vectors  are  assumed  to  be  zero  mean, 
Gaussian,  independent  random  vectors  with  identity 
covariance  (7:11). 

With  this  system  model,  the  likelihood  functional 
p(z|f)  can  be  written  directly.  However,  because  of  its 
simplicity,  the  logarithm  of  the  likelihood  ratio  is  usually 
used.  For  this  model,  the  log  likelihood  function  is 

logp(z|C)  =  -  — 

-  ^log  |GG*|  -  log  2n  Oj 


where  Zj.  is  a  predicted  estimate  of  z  computed  using  a 
Kalman  filter.  The  last  term  in  equation  (3)  is  a  constant, 
and  if  G  is  known,  the  next  to  last  term  is  also  constant. 
Therefore  they  can  be  neglected  in  the  maximization  process 
(7:11-12) . 

If  the  log  likelihood  function  is  multiplied  by  -1,  the 
problem  becomes  a  minimization  rather  than  a  maximization. 
The  resulting  problem  is  to  minimize  the  cost  function 
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Cost  Function  Minimization 


Now  that  the  cost  function  has  been  defined,  the 
algorithm  for  the  minimization  needs  to  be  specified.  The 
Newton-Raphson  technique  is  well  known  for  iterative 
functional  minimization.  The  Newton-Raphson  algorithm 
requires  an  initial  starting  point  for  the  estimated 
parameter  vector,  f.  Subsequent  values  for  f  are  determined 
from  the  equation 

Ci.i  “Ci <5) 

The  first  and  second  gradients  in  equation  (5)  are  defined 
as 

VjJ-(C)  = 

I»1 


and 


Vfj-(C)  = 

i-l 

*  ,7) 

i-l 


The  second  term  in  equation  (7)  involves  significant 
computation  for  a  realistic  number  of  parameters.  It  should 
be  small  near  the  minimum,  because  z  -  z  should  be  small. 
Therefore,  this  term  can  be  neglected  if  the  process  starts 
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close  enough  t.o  the  solution  (7:14).  This  results  in  the 
second  gradient  being  defined  by 


N 


^<7(0  = 


(8) 


This  approximation  is  called  the  Gauss-Newton 
approximation,  and  is  also  known  as  the  modified  Newton- 
Raphson  or  Newton-Balakrishnan  technique.  The  Gauss-Newton 
approximation  is  much  easier  for  computation,  since  it  does 
not  require  calculation  of  the  second  gradient  of  the 
innovations.  It  can  also  expedite  the  convergence  process 
(4:3)  . 

A-Priori  Information 

In  addition  to  the  data  available  from  processing  one 
particular  maneuver  through  the  parameter  estimation 
routine,  other  estimates  may  be  available  from  prior  runs  at 
the  same  test  case  or  from  entirely  different  sources.  The 
algorithm  can  be  set  up  to  account  for  this  information  by 
adding  a  penalty  function  to  the  cost  function  for 
deviations  from  the  a-priori  estimates.  The  quadratic 
penalty  function,  when  added  to  the  cost  function  yields  a 
cost  function  of  the  form 
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where  the  W  matrix  is  a  weighting  matrix  corresponding  to 
the  a-prior’  values  in  the  vector  fo  (7:15). 

While  the  a-priori  function  can  assist  convergence  in 
some  cases,  it  should  be  used  with  caution  since  it  also 
biases  the  estimates  (possibly  in  the  wrong  direction) .  The 
a-priori  weighting  can  be  used  only  for  initial  convergence, 
then  removed  to  arrive  at  the  unbiased  parameter  estimates. 


Cram6r-Rao  Bounds 


Since  the  maximum  likelihood  estimation  method  yields 
just  that  ...  estimates,  some  accuracy  measure  is  necessary. 
The  Cram^r-Rao  bounds  are  a  good  analytical  measure  of 
estimation  accuracy  (4:11).  In  many  cases,  the  scatter  of 
the  estimates  would  be  sufficient  for  determining  accuracy. 
However,  if  there  are  not  enough  estimates  available  for  a 
particular  test  point,  even  the  scatter  cannot  be  used.  The 
Cramer-Rao  bound  is  defined  by 


Variance^)  i  (^{[Vjlogp(z|C)]*[V(logp(z|C)]})’^ 


(10) 


It  is  important  to  note  that  the  Cramer-Rao  bounds  are  lower 
bounds  for  the  variance.  In  other  words,  the  variance  is  at 
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least  as  large  as  the  Cram^r-Rao  bound  on  the  estimate 
(7:16) . 

For  the  case  where  there  is  no  state  noise,  the  Cramfer- 
Rao  bound  can  be  determined  by 


Variance  (C) 


-1 


(11) 


Equation  (11)  is  the  inverse  of  the  Gauss-Newton 
approximation  to  the  second  gradient  from  equation  (8) . 
Because  the  second  gradient  calculation  is  already  required, 
the  Cram6r-Rao  bound  calculation  is  minimal. 

Residual  Power  Estimation 

In  most  of  the  above  discussion,  it  has  been  assumed 
that  the  G  matrix  is  known.  However,  the  G  matrix  can  be 
estimated  by  maximizing  the  log  likelihood  function  in 
equation  (3)  .  The  problem  is  that  Zj.  is  a  function  of  G. 
There  is  a  process  used  to  estimate  G  which  estimates  it  and 
the  rest  of  the  unknown  parameters  separately.  Equation  (3) 
is  maximized  with  respect  to  G,  without  considering  the 
effect  of  G  on  Zj..  In  particular,  the  estimator  is 

,12, 
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The  entire  process  Is  accomplished  by  first  estimating  the 
unknown  parameters  other  than  G.  The  G  matrix  is  then 
estimated  using  evaluated  at  the  new  ^  to  get  the  new 
estimate  of  G.  These  steps  are  repeated  until  convergence 
(7:17)  . 

Summary 

The  above  section  described  the  theory  and  algorithms 
behind  maximum  likelihood  estimation.  The  State  Space 
Identification  Tool  for  the  PC  Matlab  personal  computer 
program  implements  the  maximum  likelihood  algorithm.  The 
task  at  hand  is  to  apply  the  general  theory  to  the  specific 
application  of  aircraft  equations  of  motion.  The  process  of 
developing  models  and  getting  the  software  to  work  for 
aircraft  applications  using  both  simulated  and  flight  test 
data  is  presented  in  subsequent  sections. 
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Ill,  simulation 


Introductory  work  was  done  using  simulated  data  to 
check  out  the  personal  computer  software  for  identifying 
aircraft  response  parameters.  The  primary  areas  of  concern 
were:  developing  an  adequate  model  for  the  short  period 
aircraft  response,  generating  test  inputs  and  outputs, 
creating  the  macro  files  in  PC  Matlab  to  tailor  the  MMLE 
software  to  the  aircraft  model,  and  running  the  software 
with  the  simulated  data  to  verify  the  routines. 

Model  Development 

Two  models  were  developed  for  the  short  period  mode  of 
motion.  First,  the  short  period  approximation  was  used  to 
produce  a  differential  equation  representation  of  the  short 
period  dynamics.  The  differential  equations  were  then 
transferred  into  state  space  format  for  use  with  PC  Matlab. 
The  second  model  started  with  the  general  form  of  a  short 
period  pitch  rate  to  elevator  deflection  (q/5e)  transfer 
function.  The  primary  difference  was  that  the  transfer 
function  model  estimated  second  order  system  parameters 
(damping  ratio,  natural  frequency,  etc.)  directly,  whereas 
the  differential  equation  model  estimated  stability 
parameters  in  the  aircraft  equations  of  motion. 

Short  Period  Approximation.  The  longitudinal  small 
perturbation  equations  of  motion  can  be  expressed  in  many 
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forms.  For  simplicity  of  expression,  Roskam  (14)  develops 
the  equations  in  the  stability  axis  system  (14:391).  The 
equations  are  developed  assuming  perturbations  around  a 
steady  state  condition  for  which  there  is  no  initial  side 
velocity  or  bank  angle,  and  there  are  no  initial  angular 
velocities  (14:45).  For  convenience  the  elements  of  the 
equations  are  also  defined  in  terms  of  dimensional  stability 
parameters  (14:413).  These  equations  are  expressed  in 
differential  equation  form  (14:414)  as 

u  =  -  g0cos0i  +  X^u  +  X^u  +  X^a  + 

^  -  U^q  =  -  p0sin0j^  +  ZyU  +  Z^a  +  Z^&  +  Z^q  +  Z^^be 

and 


=  Af„u  +  Mj^u  +  W,a  +  Afy^a  +  +  M^q  + 


(15) 


For  the  short  period  approximation,  several  assumptions 
are  made  to  simplify  the  model.  First  of  all,  the  short 
period  mode  is  assumed  to  occur  at  constant  forward  speed 
(u) ,  so  the  u  portions  of  the  equations  can  be  deleted. 
Substituting  d  =  w/Uj,  the  resulting  equations  are 

{Uj^-Z^)&  =  -  grsin0i  *  (C/i+Z^)g  +  Z^^be 

and 
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<3r  =  Af,o  +  Mj.^a  +  +  M^q  +  M^^be 


(17) 


If  the  steady  state  condition  is  defined  as  trimmed 
level  flight,  then  =  0.  In  addition,  it  is  further 
assumed  that  and  Z^  are  much  smaller  than  U, .  Finally,  it 
is  assumed  that  the  derivative  is  included  in  M„.  In 
other  words,  the  combined  effect  of  and  is  included  in 
the  single  derivative  M„.  Incorporating  these  assumptions 
yields  the  following  set  of  equations: 


=  Z^a  +  U^q  +  2,^6© 


(18) 


and 


4  =  W.a  +  +  M  q  +  Afjp5e 


(19) 


One  additional  assumption  is  made  at  this  point.  While 
an  analytical  distinction  can  be  made  between  the  d  and  q 
derivatives,  flight  test  produces  poor  results  when 
attempting  to  separate  these  derivatives  in  a  particular 
maneuver  (11:17).  Therefore,  the  d  and  q  derivatives  are 
combined  into  a  single  q  derivative.  In  other  words. 


M.q  -  M.q 


(20) 
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with  this  final  assumption,  the  equations  can  now  be 
rewritten  solving  for  d  and  q.  This  results  in  the 
following  set  of  equations: 

d  =  {ZjUj)a  +  <7  +  (Z^JUj)6e 

and 


4  =  Af,a  + 


(22) 


Or,  in  state  space  representation 


11  = 

i' 

|H 

1 

>  +  < 

w 

.<7. 

1  J 

In  this  case  the  outputs  are  the  states  themselves  (angle  of 
attack  and  pitch  rate) ,  so 


y  = 


(24) 


For  the  purposes  of  checking  out  the  simulation, 
dimensional  stability  derivatives  were  used  for  the  A-4D 
aircraft.  A  middle  of  the  envelope  flight  condition  was 
chosen  at  15,000  feet  pressure  altitude  and  0.60  Mach 
number.  Table  1  below  shows  the  A-4D  data  used  (9:702). 
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Table  1.  A-4D  Siaulation  Data 


h  =  15,000  ft 

Mach  =  0.60 

Ui  =  634  ft/sec 

Z„  =  -518.6  (ft/sec^) /rad 

M„  =  -12.68  (l/sec^) 

Z^e  =  -56.92  (ft/sec^) /rad 

M,,  =  -1.07  (1/sec) 

Mfc  =  -19.48  (1/sec^) 

Rewriting  equation  (23)  with  the  specific  value  of  U,  yields 


(*1 

(.00158 

2.)  l' 

|«1 

^  |(.  00158  2,4 

^1 

t6e 


(25) 


Now  for  purposes  of  parameterizing  tne  model  for  the 
MMLE  subroutine,  the  stability  derivatives  in  equation  (25) 
are  represented  as  follows: 

P(l)  =  Z„ 

P(2)  = 

P(3)  =  M, 

P(4)  = 
p(5)  =  Mfc 

So,  the  final  form  of  the  model  for  the  MMLE  subroutine  is 


H  - 

.00158*p(l) 

1 

f.00158*p(4) 

P(2) 

P(3)  , 

m 

1  P(5) 
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y  = 


1  0 
0  1 


+ 


(27) 


Equations  (26)  and  (27)  represent  the  final  form  of  the 
short  period  approximation  model  which  was  used  for  the  MMLE 
software  checks  using  simulated  data. 

Transfer  Function  Model.  In  addition  to  the  short 
period  approximation  model,  a  second  method  was  used  to 
provide  a  cross  check  of  the  results  from  the  short  period 
approximation.  This  second  method  estimated  second  order 
system  parameters  directly,  and  the  model  development 
started  with  the  short  period  approximation  pitch  rate  to 
elevator  deflection  (q/5e)  transfer  function: 


g{s) 
be  is) 


s2+2C5pa>^pS+co|p 


(28) 


A  set  of  state  variables  is  required  so  that  this  transfer 
function  can  be  expressed  in  the  form  x  =  Ax  +  Bu.  Since 
the  transfer  function  could  be  rewritten  in  differential 
equation  form,  the  resulting  differential  equation  was 
examined  to  help  define  state  variables.  Equation  (28)  can 
be  rewritten  as 


19 


(s2  +  2Cgp«spS  +  <olp)g(s)  =  Ke{s  +  l/re2)6e(s) 


(29) 


Taking  the  inverse  Laplace  transform  yields 

g  +  +  «|pg  =  K^6e  +  K^{UT^2)^e  (30) 

Equation  (30)  is  the  differential  equation 
representation  of  the  short  period  q/^e  transfer  function. 
Normally  q  and  q  could  be  used  for  state  variables  since 
knowledge  of  q(0)  and  q(0)  along  with  5e(t)  for  t  >  0 
completely  determines  the  response  of  the  system.  However, 
a  problem  arises  due  to  the  presence  of  the  derivative  term 
of  the  forcing  function  5e(t)  in  equation  (30) .  Ogata  (12) 
presents  a  method  of  choosing  state  variables  whereby  the 
control  derivative  term  is  eliminated  from  the  state 
equation  (12:675-678). 

If  a  set  of  state  variables 


and 


X,  =  q 

X2  =  q  -  K(,5e 


is  chosen,  then  the  state  and  output  equations  for  the 
differential  equation  in  equation  (30)  are 
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=  4r  =  Xj  + 


(31) 


^2  =  ^  -  2Csp“sp'3r  +  k»{^/Tq2)^^ 


=  -olpXi  -  2C^pa)^pX2  +  A^(l /  Tea  -  2  CspWsp)  5  e 


(32) 


and 


y  =  Xi 


(33) 


Rewritten  in  matrix  form,  equations  (31)  through  (33)  become 


0  1 

k 

1 J  ^  1 

-2C  (•) 

«P  “  » flp  sp 

\x, 

1  K(l/rea-2C.p0.p)J 

(34) 


and 


y  = 


[1  0] 


(35) 


Now  to  parameterize  the  model  for  MMLE,  the  terms  in 
equation  (34)  are  redefined  as  follows; 

P(l)  =  K<, 
p(2)  =  1/T„ 
p(3)  = 
p(4)  = 
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Therefore,  the  final  model  used  for  the  pitch  rate  transfer 
function  model  with  MMLE  is 


0 

-p(4) 


1 

-P(3) 


+ 


I  P(l) 

[p(l)  *  {p(2)  -p(3)} 


(36) 


and 


y  =  g  = 


[1  0] 


(37) 


For  the  purposes  of  checlcing  out  the  simulation,  the 
transfer  function  parameters  were  calculated  based  on  the 
dimensional  stability  derivatives  for  the  A-4D  aircraft. 
These  derivatives  were  shown  in  Table  1  above.  Starting 
with  the  Laplace  transforms  of  equations  (21)  and  (22) , 

sa(s)  =  (Z^/Uj)a(s)  +  q(s)  +  (ZiJUj)6e{s)  ^3 

and 


sgis)  =  M^a(s)  +  Mgq(s)  +  M^^beis) 


(39) 


or  in  matrix  form. 


'{3-Za/U,) 

-1 

a  (s) 

-M, 

{s-M^y 

q(s) 

.  ) 


be(s) 


(40) 
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Using  Cramer's  rule  on  equation  (40)  yields 


g(g)  ^  *  M^jZ.JU,) 

5e(s)  [s-{ZjUj)]{s-Mg)  -  M,  (41) 


Equation  (41)  can  be  rewritten  as 

g(s)  =  *  \M^(Z,JU^)-Mi,(ZjU^)] 

5e(s)  s2  -  [(Z,/Uj)*M^]s  +  [M^Z,/U^)-M^]  (42) 

Substituting  the  A-4D  data  from  Table  1, 

q(s)  ^  -19.48  is  +  0.75954) 

6e(s)  +  1.88798s  +  13.55524  (43) 


Equation  (43)  is  now  in  the  standard  form  of  equation  (28) . 
The  transfer  function  parameters  of  interest  (Kj,  l/T^j, 
and  «^)  can  be  calculated  directly  from  equation  (43)  .  The 
parameters  calculated  using  the  A-4D  data  are  shown  in  Table 
2  below.  These  numbers  will  be  used  to  compare  with  the 
MMLE  estimates  obtained  using  the  transfer  function  model. 


Table  2.  A-4D  Calculated  Transfer  Function  Data 


h  =  15,000  ft 

Mach  =  0.60 

00 

• 

0\ 

1 

II 

1/T„  =  0.7595 

=  0.2564 

=  3.6817 
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Input  and  Output  Generation 


The  type  of  input  used  for  the  parameter  estimation 
maneuver  is  critical  to  the  success  of  the  process.  Past 
studies  have  shown  that  inputs  with  sharp  corners  produced 
the  best  results.  Sets  of  derivatives  have  been  obtained 
using  simple  doublet  inputs,  which  produce  results  as  good 
as  those  from  more  complex  inputs  (15:1). 

For  the  longitudinal  axis,  the  maneuver  of  choice  at 
the  AFFTC  is  a  rapid  elevator  doublet.  The  doublet  is 
followed  by  approximately  5  to  10  seconds  of  free  response. 
The  control  derivatives  are  identified  during  the  input  time 
segment,  and  stability  and  damping  derivatives  are 
identified  during  the  free  response. 

As  described  by  Nagy  (11:31),  derivatives  in  general 
are  functions  of  angle  of  attack  and  Mach  number,  and  these 
are  the  key  parameters  which  should  be  held  as  constant  as 
possible  throughout  the  maneuver.  Although  angle  of  attack 
will  obviously  vary  with  an  elevator  doublet,  the  doublet 
should  be  short  enough  in  duration  that  no  large  pitch  rate 
or  angle  of  attack  perturbations  occur.  Then  at  least  the 
derivatives  will  be  valid  for  the  average  angle  of  attack  of 
the  maneuver.  As  for  Mach  number,  the  importance  of 
maintaining  trim  conditions  is  a  primary  reason  for  using  a 
doublet  input  instead  of  a  pulse  (11:31). 
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The  doublet  input  should  have  a  magnitude  such  that  the 
response  is  as  large  as  possible  without  violating  the 
assumptions  of  the  linear  model.  Larger  inputs  result  in 
higher  signal-to-noise  ratios  in  the  data,  and  hence  yield 
better  results  (15:6).  As  a  rough  guideline,  pitch  rate 
excursions  should  be  around  10  to  20  degrees  per  second  peak 
to  peak.  Angle  of  attack  should  vary  within  about  four 
degrees  peak  to  peak  (11:31). 

A  triangular  doublet  was  designed  using  PC  Mat lab,  with 
variable  amplitude  and  frequency.  For  this  simulation,  the 
frequency  was  chosen  to  give  a  complete  doublet  input  in  one 
second.  The  amplitude  was  chosen  as  ±0.07  radians  of 
elevator  deflection  (approximately  ±4  degrees) .  The  input 
was  designed  with  20  samples  per  second  to  avoid  problems  in 
determining  the  effects  of  an  input  at  this  frequency.  The 
total  maneuver  duration  was  10  seconds. 

The  "true”  values  of  each  parameter  (shown  in  Table  1) 
for  the  A-4D  were  used  to  generate  the  output  of  the  system 
With  the  Table  1  values,  the  resulting  state  space  system  is 


-.81939  1  |al  f-.08993K 

-12.68  -1. 07  Jlgj  ^  1-19.48  j  ® 


(44) 


y  = 


(45) 
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The  noise  free  output  response  data  for  pitch  rate  and 
angle  of  attack  are  generated  using  the  DLSIN  function  in  PC 
Mat lab  (6:CR-32) .  The  DLSIM  function  simulates  the  outputs 
of  a  discrete  time  linear  system.  The  state  space  system  in 
equations  (44)  and  (45)  was  discretized  using  the  C2D 
continuous  to  discrete  conversion  function  of  PC  Matlab 
(6:CR-11).  The  C2D  function  requires  the  system  A  and  B 
matrices,  and  the  time  increment  as  inputs.  It  returns  the  i 
and  r  matrices  in  the  corresponding  discrete  time  system: 

x[n+l]  =  <tx[n]  +  ru[n] 

assuming  the  control  input  is  piecewise  constant  over  the 
time  increment.  The  DLSIM  function  takes  the  ♦  and  r 
matrices,  along  with  the  output  C  and  D  matrices  and  the 
defined  input,  and  returns  the  simulated  time  history  data 
for  th3  specified  system  and  input. 

The  state  space  system  matrices  of  equations  (44)  and 
(45)  were  processed  through  the  continuous  to  discrete 
conversion  function,  and  subsequently  through  the  discrete 
linear  systems  simulation.  Figure  1  shows  the  PC  Matlab 
code  required  to  define  the  input,  input  the  system 
matrices,  and  produce  the  time  history  output.  Figure  2 
shows  the  simulated  time  histories  of  elevator  deflection, 
angle  of  attack,  and  pitch  rate,  respectively.  The  noise 
free  simulated  data  was  then  used  to  accomplish  the 
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disp(ans) 

dispC  INPUT/OUTPUT  SIMULATION  MACRO*) 

dispC  A-4D  LONGITUDINAL  STABILITY  PARAMETERS') 

dispC  TRUE  DATA  =  [-518.6  -12.68  -1.07  -56.92  -19.48]') 
disp(ans) 

% 

% - SIMULATE  INPUT 

% 

ndp=201;  %  INPUT  NUMBER  OF  DATA  POINTS 

dt=.05;  %  INPUT  TIME  INCREMENT 

ainp=.07;  %  INPUT  AMPLITUDE  OF  DOUBLET  IN  RADIANS 

per=l;  %  INPUT  PERIOD  OF  DOUBLET  IN  SECONDS 


ndp=201;  %  INPUT  N1 

dt=.05;  %  INPUT  t: 

ainp=.07;  %  INPUT  A] 

per=l;  %  INPUT  PI 

% 

t= [ 0 : ndp-1 ] *dt ; 
uydata=zeros (ndp, 3) ; 
inc= (4*amp*dt) / (per) ; 
ul=(0:-l*inc;-l*ainp)  ; 
u2=  ( (-l*ainp)  +inc:  inc:  amp)  ; 


%  THESE  LINES 
%  GENERATE  THE 
%  ELEVATOR  DOUBLET 
%  FROM  THE  USER'S 
%  SPECIFIED  PERIOD 
%  AND  AMPLITUDE 


u3=(amp-inc:-l*inc:0) ;  %  AND  AMPLITUDE 

uydata ( : , 1) =[ul  u2  u3  zeros (1, ndp- (per/dt) -1) ]' ; 

% 

% - SIMULATE  OUTPUT 

% 

%  INPUT  SYSTEM  MATRICES  FOR  SIMULATION 

% 

a=[-. 81939  1 

-12.68  -1.07]; 

% 

b=[-. 08993 
-19.48  ]; 

% 


UYDATA (:, 2)  =  ANGLE  OF  ATTACK  (RAD) 
UYDATA(;,3)  =  PITCH  RATE  (RAD/SEC) 


[phi,gam]=c2d(a,b,dt) ;  %  DISCRETIZE 
% 

uydata  ( : ,  2 : 3)  =dlsim  (phi,  gam,c,d,  uydata  { : ,  1) )  ;  %  NOISE-FREE 
save  uydata  %  RESPONSE 
% - ENDSIMULATE.M 


Figure  1.  PC  Natlab  Data  Simulation  Code 
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Flgur«  2.  simulated  Data  Time  Histories 
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parameter  estimation  task  through  both  the  short  period 
approximation  and  transfer  function  models. 

PC  Matlab  Macro  File  Creation 

Using  the  PC  Matlab  MMLE  program  to  accomplish 
parameter  identification  involves  basically  three  steps. 
First,  the  user  must  create  a  file  to  convert  the  parameter 
vector  into  a  state  space  system  model.  Secondly,  the 
necessary  inputs  to  the  MMLE  subroutine  must  be  defined. 
Finally,  the  MMLE  macro  file  is  executed. 

The  first  step  is  the  creation  of  the  parameter-to- 
state-space  (P2SS)  macro  function.  This  function  primarily 
converts  the  parameter  vector  (p)  into  the  user-defined 
state  space  description  of  the  system.  In  addition  to  the 
state  space  system  matrices  (A,  B,  C,  and  D) ,  the  function 
must  return  the  following; 

1.  q  =  square  root  of  the  state  noise  covariance 
matrix. 

2.  rowing  =  a  vector  (the  same  size  as  p)  showing 
where  parameters  to  be  estimated  occur  in  the  q 
matrix. 

3 .  xO  =  a  row  vector  containing  initial  conditions  of 
the  states. 

4.  dt  =  the  time  increment  for  discretizing  the 
continuous  model. 
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5.  phi  =  9,  the  state  transition  matrix. 

6.  gam  =  r,  the  discrete  control  distribution  matrix. 

The  i  and  r  matrices  can  be  obtained  from  the  continuous  to 
discrete  transformation  (C2D)  already  discussed. 

Once  the  P2SS  file  is  created,  several  other  necessary 
inputs  to  MMLE  must  be  defined.  The  easiest  way  is  to 
create  another  macro  file  for  defining  the  required  inputs 
to  MMLE.  The  following  are  the  necessary  inputs  (10:SR- 
6,SR-7) : 

1.  uydata  =  a  matrix  consisting  of  columns 
corresponding  to  the  time  histories  of  first  the 
inputs,  then  the  outputs. 

2.  p2snam  =  a  string  containing  the  file  name  of  the 
P2SS  function. 

3.  pO  =  a  row  vector  containing  initial  values  of  the 
parameters  to  be  identified. 

4.  rmso  =  a  vector  (the  same  size  as  pO)  containing 
standard  deviations  of  a-priori  estimates  if  the 
a-priori  weighting  option  is  used.  The  rmsO  value 
defaults  to  zero  if  it  is  undefined. 

5.  pref  =  a  vector  (the  same  size  as  p)  containing 
reference,  or  true,  values  for  each  of  the 
parameters . 
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6.  pidq  =  a  row  vector  giving  the  locations  of 
parameters  in  the  vector  p  which  should  be 
identified  in  the  quadratic  phase. 

7.  pidm  =  a  row  vector  giving  the  locations  of 
parameters  in  the  vector  p  which  should  be 
identified  in  the  Marguardt  phase. 

8.  pidf  =  a  row  vector  giving  the  locations  of 
parameters  in  the  vector  p  which  should  be 
identified  in  the  constrained  Newton  and  "gg 
determination"  phases. 

9.  pert  =  the  perturbation  to  be  used  for  numerical 
calculation  of  gradients. 

10.  ggO  =  the  initial  guess  of  the  innovations 
covariance  matrix. 

11.  diagrr  =  a  flag  which,  if  defined,  zeroes  off- 
diagonal  terms  in  the  estimated  innovation 
covariance  matrix. 

12.  linesearch  =  1  enables  the  linesearch  along  the 
Newton  direction.  Linesearch  defaults  to  zero  if 
it  is  not  defined. 

13.  opt  =  a  vector  which  controls  the  number  of 
iterations  and  criteria  for  convergence.  The  opt 
default  is  [0  5  5  5  .02  .05  .001  1].  The 
elements  of  the  vector  (in  sequence)  are  as 
follows: 
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a.  qstep  =  the  number  of  quadratic  steps. 

b.  mstep  =  maximum  number  of  Marquardt  steps. 

c.  nstep  =  maximum  number  of  constrained  Newton 
steps. 

d.  gstep  =  maximum  number  of  "gg  determination” 
steps . 

e.  marqf  =  the  value  of  parameter  marq  which 
terminates  the  Marquardt  phase. 

f.  fconv  =  the  convergence  criterion  for  early 
exit  from  the  constrained  Newton  phase. 

g.  gconv  =  final  convergence  criterion. 

h.  rhz  =  frequency  of  a  low  pass  filter  used  to 
compute  filtered  Cramer-Rao  bounds. 

The  four  minimization  phases  discussed  above  require 
further  description  (10:SR-7,SR-8) .  The  quadratic  phase  is 
a  pure  Newton  phase  which  can  find  the  minimum  of  a 
quadratic  cost  function  in  one  iteration.  This  phase  is 
skipped  in  the  default  opt  vector.  The  Marquardt  phase 
varies  the  step  size  along  the  gradient  to  ensure  that  the 
problem  is  well  conditioned  for  the  constrained  Newton 
phase.  The  constrained  Newton  phase  is  the  final 
optimization  phase  which  uses  the  Newton  algorithm.  The  ”gg 
determination”  phase  is  used  for  estimating  the  gg  matrix. 
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Once  the  necessary  Inputs  are  defined,  completing  the 
parameter  identification  process  is  simply  a  matter  of 
executing  the  MMI^S  macro.  Having  addressed  the  requirements 
for  running  MMLE  in  PC  Matlab,  the  specific  macro  files  for 
accomplishing  parameter  identification  with  the  two  models 
can  be  discussed. 

Short  Period  Approximation  Macro.  Recall  from 
equations  (26)  and  (27)  the  final  form  of  the  short  period 
approximation  model  parameterized  for  use  with  MMLE.  The 
equations  are  repeated  here  for  convenience: 


H  = 

’.00158*p(l)  1 

/a1  /.00158>*p(4)' 

laj 

P(2)  p(3). 

W  1  P(5) 

(47) 


and 


y  = 


1  0 
0  1 


a 


?  + 


(48) 


The  macro  file  SPP2SS.M  was  created  to  perform  the 
parameter-to-state-space  conversion  for  the  short  period 
approximation  model.  The  PC  Matlab  code  in  the  macro  file 
is  shown  in  Figure  3.  The  matrix  definitions  are  straight 
forward.  Since  the  MMLE  run  is  to  be  performed  against 
simulated  data,  the  state  noise  matrix  Q  is  set  to  a  small 
non-zero  diagonal  matrix  (basically  eliminating  its  effect 


33 


f unct i on [ a , ph i , gam , c , d , q , xO , dt , rowing , b ] =spp2  s s ( p ) ; 

% - 

%  P2SS  FUNCTION  FOR  SHORTPD.M  -  EXTENSIVELY  COMMENTED 

%  STATE-SPACE  REALIZATION  OF  SHORT  PERIOD  APPROXIMATION 
%  WITH  A-4D  STABILITY  PARAMETERS  (H  =  15,000  FT,  M  =  0.6) 

% - 

% 

%  p(l)  =  Z_ALPHA 
%  p(2)  =  M_ALPHA 


%  p(3)  =  M_Q 
%  p(4)  =  Z_DE 
%  p(5)  =  M_DE 
% 

% - SYSTEM  MATRICES 

% 

a=[ .00158*p(l)  1 

P(2)  P(3)]; 

% 

b=t .00158*p(4) 


P(5)  ]; 

% 

c=[  1  0 

0  1]; 

% 

d=[  0 

0  ]; 

% 


% - STATE  NOISE  COVARIANCE 

% 

q=eye(2) *le-6;%  Q  IS  THE  SAME  SIZE  AS  A  WITH 

%  Q*Q'  POSITIVE  DEFINITE 1 

% 

% - ROWS  IN  Q  IN  WHICH  PARAMETERS  OCCUR 

% 

rowinq=0*p; 

% 

% - INITIAL  STATE  VECTOR 

% 

x0=zeros(l,2) ; 

% 

% -  DEFINE  DT  AND  DISCRETIZE 

% 

dt=. 05; 


[phi,gam]=c2d(a,b,dt) ; 

% 

% - end  SPP2SS.M 

Figure  3.  short  Period  P2SS  Function  Code 
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on  the  problem) .  There  are,  therefore,  no  parameters  to  be 
identified  in  the  Q  matrix,  so  the  rowing  variable  is  a 
vector  of  zeroes.  The  simulated  data  started  with  zero 
initial  conditions,  so  the  xO  vector  is  also  set  to  zero. 
Finally,  the  macro  file  defines  the  time  increment  from  the 
data  simulation  and  makes  the  call  to  the  continuous  to 
discrete  (C2D)  function  to  calculate  the  *  and  F  matrices. 

Once  the  SPP2SS.M  file  is  created,  the  required  inputs 
discussed  above  have  to  be  defined.  This  purpose  is  served 
by  the  SHORTPD.M  file,  shown  in  Figure  4.  This  macro  is  set 
up  to  input  the  true  parameters  for  the  model  (previously 
shown  in  Table  1) ,  and  load  the  input  and  output  data  file 
created  by  the  simulation  macro  defined  above.  The  pidq, 
pidffl,  and  pidf  parameters  are  defined  such  that  all  five 
parameters  for  the  identification  will  be  identified  in  the 
constrained  Newton  phase.  Even  though  the  quadratic  and 
Marquardt  phases  are  not  used,  the  corresponding  parameters 
(pidq  and  pidffl)  must  be  defined.  Therefore,  they  are  set  to 
a  nominal  value  of  one,  while  pidf  is  set  to  reflect 
identification  of  parameters  one  through  five. 

The  initial  parameter  estimate  is  set  up  with  a  30 
percent  error  from  the  true  parameters.  The  default  opt 
vector  is  redefined  to  reflect  parameter  identification  in 
the  constrained  Newton  phase  only,  and  has  a  slightly  more 
stringent  requirement  on  final  convergence.  The  ggO  matrix 
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format  compact, clc 
clear; 

diary  shortpd.log 


I  _ _ _ I . 

dispC  SHORT  PERIOD  APPROXIHATION  MODEL') 

dispC  A-4D  LONGITUDINAL  STABILITY  PARAMETERS') 

dispC  TRUE  DATA  =  [-518.6  -12.68  -1.07  -56.92  -19.48]') 
% 

% - GET  TRUE  STABILITY  PARAMETERS 

% 


%  THE  RESULTS  MATRIX  PRODUCED  BY  MMLE.M  USES  PREF  IF 
%  AVAILABLE 
% 

pref=[-518.6  -12.68  -1.07  -56.92  -19.48]; 


% 

% - SUFFICIENT  INPUTS  TO  MMLE  FOLLOW 

% 

load  uydata;  %  UYDATA  IS  OUTPUT  FROM  THE  SIMULATE. M  MACRO, 
%  AND  MUST  BE  AVAILABLE  AS  A  DATA  FILE. 

% 

pidq=[l];  %  IDENTIFY  WHICH  PARAMETERS  ARE  TO  BE  ESTIMATED 
pidm=[l];  %  IN  THE  QUADRATIC  (PIDQ) ,  MARQUARDT  (PIDM) , 


pidf=[l:5];%  AND  CONSTRAINED  NEWTON  (PIDF)  PHASES 
% 

p0=0.7*pref  ;%OFFSET  INITIAL  GUESS  30%  FROM  TRUE  PARAMETERS 

% 

opt=[0  0  5  0  .02  .01  .001  1];  %IF  YOU  SIMULATED  NOISE-FREE 
%  DATA,  MAKE  SURE  OPT(4)=0! 

% 

gg0=le-4*eye(2) ;  %  INITIAL  GUESS  OF  INNOVATIONS  COVARIANCE 
% 

p2snam='spp2ss' ;  %  FILENAME  OF  P2SS  FUNCTION 

% 

pert=le-4;  %  SET  PERTURBATION  FOR  GRADIENT  CALCULATION 
% 

linesearch=l ; 

mmle  %  SUBROUTINE  CALL  TO  MMLE  PROGRAM 

% 

Z_ALPHA=pfin(l)  %  PRINT  OUT  FINAL  STABILITY  PARAMETERS 
M_ALPHA=pfin(2) 

M_Q=pfin(3) 

Z_DE=pfin(4) 

M_DE=pfin(5) 

% 

diary  off 

% - end  shortpd.m 

Figure  4.  Short  Period  MMLE  Macro  Code 
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is  defined,  as  with  Q,  as  a  small  non-zero  valued  diagonal 
matrix.  Finally,  the  linesearch  flag  is  set  and  the 
subroutine  call  to  MMLE  is  made.  The  remaining  lines  in  the 
macro  file  are  for  final  printout  of  the  parameter 
identification  results. 

Transfer  Function  Macro.  The  final  form  of  the  transfer 
function  model  for  use  with  the  MMLE  program  was  expressed 
in  equations  (36)  and  (37) ,  and  is  repeated  here  for 
convenience: 


0  1 
-p(4)  -p(3) 


H  r  , 

Ipd)  *  {p(2)  -p(3)}  J 


(49) 


and 


y  =  [1  0] 


(50) 


For  the  transfer  function  model,  the  QDEP2SS.M  file 
sets  up  the  parameter-to-state-space  conversion.  As  with 
the  short  period  approximation  model,  the  file  defines  the 
state  space  matrices  in  terms  of  the  parameter  vector,  and 
the  remaining  variable  definition  is  identical  to  that  in 
the  SPP2SS.M  file.  The  QDEP2SS.M  macro  file  is  shown  in 
Figure  5. 

As  with  the  short  period  approximation  model,  a  macro 
file  performs  all  of  the  required  input  definition  for  MMLE. 
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f unction [a, phi, gam, c,d,q,xO,dt, rowing, b]=qdep2ss(p) ; 

% - 

%  P2SS  FUNCTION  FOR  QDELOES.M  -  EXTENSIVELY  COMMENTED 

%  STATE-SPACE  PITCH  RATE  TO  ELEVATOR  X-FER  FCN 


% - 

%  Q  A_THETA  *  (S  +  1/T_THETA2) 

%  ___  =  - 

%  dE  S'‘2  +  2*ZETA*WN*S  +  WN^2 

% - 

% 

%  p(l)  =  A_THETA 
%  p(2)  =  1/T_THETA2 
%  p(3)  =  2*ZETA*WN 
%  p(4)  =  WN^2 
% 

% - SYSTEM  MATRICES 

% 

a=[  0  1 

-l*p(4)  -l*p(3)]; 

% 

b=[  p(l) 

p(l)*(p(2)-p(3))]; 

% 

c=[  1  0  ]; 

% 

d=0; 

% 

% - STATE  NOISE  COVARIANCE 

% 

q=eye(2) *le-6;%  Q  IS  THE  SAME  SIZE  AS  A  WITH 

%  Q*Q'  POSITIVE  DEFINITE! 

% 

% - ROWS  IN  Q  IN  WHICH  PARAMETERS  OCCUR 

% 

rowinq=0*p; 

% 

% - INITIAL  STATE  VECTOR 

% 

x0=zeros(l,2)  ; 

% 

% - DEFINE  DT  AND  DISCRETIZE 

% 


dt=.05; 

[phi,gam]=c2d(a,b,dt) ; 

% 

% - ENDQDEP2SS.M 


Figura  5.  Transfer  Function  P2SS  Function  Code 
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The  QDELOES.M  file  (shown  in  Figure  6)  defines  all  the 
necessary  variables,  with  the  same  assumptions  and 
definition  used  in  the  SHORTPD.M  macro  file  used  for  the 
short  period  approximation  model.  The  primary  difference 
between  the  two  is  that  the  QDELOES.M  file  calculates  the 
equivalent  second-order  system  parameters  (Kj,  l/T^j/  f,p,  and 
<j}^)  ,  and  requires  only  pitch  rate  as  an  output  parameter. 

Simulation  Results 

The  simulated  input  and  output  data  described  above 
were  processed  through  both  the  short  period  approximation 
and  the  transfer  function  models.  Since  the  data  were 
simulated,  and  produced  responses  with  no  noise  content,  one 
would  expect  the  MMLE  program  to  be  able  to  match  the  time 
history  responses  exactly.  That  is,  in  fact,  what  happened. 
For  both  models,  the  predicted  time  history  data  from  the 
MMLE  program  were  identical  to  the  "actual"  simulated  time 
histories. 

Since  the  data  were  simulated  using  the  data  for  the  A- 
4D  in  the  short  period  approximation  model,  one  would  also 
expect  the  parameters  identified  to  match  exactly.  That  was 
also  the  case,  within  very  tight  tolerances.  The  A-4D 
parameters  were  listed  in  Table  1.  They  are  compared  with 
the  results  from  the  short  period  approximation  run  in  Table 
3  below. 
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format  compact, clc 
clear; 

diary  qdeloes.log 


• _ _ _ I . 

dispC  PITCH  RATE  TO  ELEVATOR  EQUIVALENT  SYSTEM  MATCH') 
dispC  A-4D  LONGITUDINAL  SYSTEM  PARAMETERS') 

dispC  TRUE  DATA  =  [-19.48  .75954  1.88798  13.55524]') 

% 

% - SIMULATE  DATA 

% 


%  THE  RESULTS  MATRIX  PRODUCED  BY  MMLE.M  USES  PREF  IF 
%  AVAILABLE 
% 

pref=[-19.48  .75954  1.88798  13.55524];  %  GET  TRUE  SYSTEM 
%  PARAMETERS 

% 

% - SUFFICIENT  INPUTS  TO  MMLE  FOLLOW 

% 

uydata=  [  uydata  ( :  1 )  uydata  (  : ,  3 )  ]  ;  %  USES  INPUT/OUTPUT  VECTOR 
%  CREATED  BY  SHORT  PERIOD 

%  SIMULATION,  WITH  ONLY 

%  PITCH  RATE  AS  AN  OUTPUT 

% 

pidq=[l];  %  IDENTIFY  WHICH  PARAMETERS  ARE  TO  BE  ESTIMATED 
pidm=[l:4];%  IN  THE  QUADRATIC  (PIDQ) ,  MARQUARDT  (PIDM) , 
pidf=[l;4];%  AND  CONSTRAINED  NEWTON  (PIDF)  PHASES 
% 

p0=0.7*pref ;%OFFSET  INITIAL  GUESS  30%  FROM  TRUE  PARAMETERS 
% 

opt=[0  0  5  0  .02  .01  .001  1];%  IF  YOU  SIMULATED  NOISE-FREE 
%  DATA,  MAKE  SURE  OPT(4)=0! 

% 

gg0=le-4;  %  INITIAL  GUESS  OF  INNOVATIONS  COVARIANCE 

% 

p2snam='qdep2ss' ;  %  FILENAME  OF  P2SS  FUNCTION 

% 

pert=le-4;  %  SET  PERTURBATION  FOR  GRADIENT  CALCULATION 
linesearch=l; 

mmle  %  SUBROUTINE  CALL  TO  MMLE  PROGRAM 

% 

a_theta=pfin(l)  %  PRINT  OUT  FINAL  SYSTEM  PARAMETERS 

invt_theta2=pf in ( 2 ) 

zetasp=pfinp)  /  (2*sqrt(pfin(4) ) ) 

wsp=sqrt ( pf in ( 4 ) ) 

diary  off 

% - ENDQDELOES.M 


Figure  6.  Transfer  Function  MMLE  Macro  Code 
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Tabl«  3.  Short  Poriod  MMLE  Rosults 


PARAMETER 

"TRUE"  DATA 

MMLE  DATA 

z„ 

-518.60 

-518.6013 

-12.68 

-12.6800 

-1.07 

-1.0700 

-56.92 

-56.9177 

1  Mfc 

-19.48 

-19.4800 

The  MMLE  run  using  the  transfer  function  model  also 
matched  the  time  histories  exactly.  The  parameters 
identified  using  this  model  were  compared  with  those 
calculated  using  Cramer's  rule  on  the  short  period 
approximation  system  matrices.  Those  calculations  were 
listed  in  Table  2.  Table  4  below  summarizes  the  MMLE 
results  from  the  transfer  function  model. 


Table  4.  Transfer  Function  MMLE  Results 


PARAMETER 

"TRUE"  DATA 

MMLE  DATA 

K, 

-19.4800 

-19.4800 

1/T« 

0.7595 

0.7609 

0.2564 

0.2566 

3.6817 

3.6819 

Having  achieved  such  good  results  with  the  simulated 
data,  the  MMLE  routines  were  considered  to  be  checked  out 
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and  ready  to  use  with  flight  test  data.  However,  several 
issues  surfaced  prior  to  flight  test  which  required 
modification  of  the  longitudinal  model.  The  next  section 
discusses  the  reasons  behind  the  changes,  and  the  subsequent 
model  which  was  developed  for  use  in  flight  test. 
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IV.  Preparations  for  Flight  Test 


Model  Changes 

Several  issues  arose  prior  to  flight  test  which  brought 
about  changes  in  the  model.  The  primary  reasons  were 
violated  assumptions  and  convenience  for  flight  test 
parameter  measurements.  The  specifics  behind  the  changes 
are  discussed  in  this  subsection. 

As  previously  mentioned,  the  equations  of  motion  used 
for  the  short  period  approximation  model  were  derived  with 
several  assumptions.  First  of  all,  they  were  developed  in 
terms  of  dimensional  stability  parameters.  Non-dimensional 
stability  derivatives  are  the  preferred  form,  because  they 
facilitate  comparisons  among  different  types  of  aircraft,  as 
well  as  different  flight  conditions.  A  form  of  the  model 
using  these  non-dimensional  stability  derivatives  directly 
was  considered  preferable  to  that  developed  for  the 
simulation. 

The  previous  model  was  also  developed  using  the 
stability  axis  system.  The  parameters  measured  in  flight 
test  are  normally  measured  with  respect  to  the  body  axis  of 
the  aircraft.  In  addition,  the  data  available  for  moments 
and  products  of  inertia  are  normally  referenced  to  the  body 
axis.  While  the  transformation  from  stability  to  body  axis 
is  a  simple  rotation  through  the  aircraft  angle  of  attack, 
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it  was  considered  preferable  to  express  the  model  in  terms 
of  the  aircraft  body  axis  to  simplify  the  calculations 
required. 

The  non-linear  equations  of  motion  were  linearized 
about  a  steady  state  condition  of  trimmed,  wings  level 
flight.  One  of  the  obvious  assumptions  associated  with  this 
condition  is  that  there  is  no  initial  bank  angle.  For 
performing  the  required  elevator  doublets  from  trimmed, 
wings  level  flight,  this  assumption  would  be  valid. 

However,  as  will  be  discussed  in  the  flight  test  section, 
one  of  the  methods  for  performing  these  maneuvers  includes 
steady  state  conditions  of  trimmed,  elevated  load  factor 
turns.  These  conditions  are  used  to  obtain  flight  test  data 
at  incremental  angles  of  attack  above  trim,  while 
maintaining  a  particular  flight  condition.  In  this  case, 
the  assumption  of  no  initial  bank  angle  is  violated,  and  the 
relationships  among  parameters  in  the  equations  of  motion 
are  more  complex. 

From  the  above  issues,  it  can  be  seen  that,  whether  out 
of  convenience  or  necessity,  changes  are  required  in  the 
model  prior  to  flight  test.  Since  the  MMLE3  program  is 
already  accepted  as  the  flight  test  standard  for  parameter 
identification,  it  is  a  good  starting  point  for  the 
longitudinal  aircraft  model. 
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The  New  Model 


The  longitudinal  equations  of  motion  developed  in  the 
body  axis  for  use  in  MMLE3  are  (11:206;7:44) 


&  =  q  -  ptanPcosa  + 


gcos0cos<|> 

VcosacosP 


qS 

WcosacosP 


Mb'- 


(51) 


-‘■yy 

+  (±-pq)  +  qSc  ( +C^^5  e+  ( C^q) )  ]  ( 52 ) 


and 


•  « 

0  =  qcos^  -  rsinij)  +  0^  (53) 

Since  the  state  equations  contain  terms  which  are  neither 
functions  of  the  state  variables  nor  the  control  input,  the 
control  vector  is  augmented  with  an  additional  input  to 
incorporate  these  terms.  It  should  also  be  pointed  out  that 
the  terms  such  as  C^o,  C^^,,  and  §„  in  these  equations  are 
also  used  to  correct  for  instrumentation  measurement  biases. 
Some  problems  were  encountered  in  the  initial  development  of 
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the  model,  when  considering  the  variables  as  perturbations 
from  trim.  If  these  terms  are  not  included  in  the  model, 
any  instrumentation  biases  translate  directly  into  initial 
rates  in  the  model.  This  problem  is  overcome  by  maintaining 
the  variables  as  total  measurements,  with  the  appropriate 
bias  terms  included  in  the  model. 

Equations  (51)  through  (53)  represent  a  two  degree  of 
freedom  model,  where  the  forward  velocity  equation  is 
omitted.  The  derivatives  in  that  equation  are  normally 
obtained  from  other  sources.  These  equations  are  non¬ 
linear,  and  several  assumptions  must  be  made  to  reduce  them 
to  a  form  usable  by  the  MMLE  routine.  It  is  first  assumed 
that  there  is  negligible  coupling  with  the  lateral- 
directional  axes.  In  other  words,  for  a  pure  elevator 
doublet,  p  =  r  =  0.  It  is  also  assumed  that  the  x-z  plane 
is  a  plane  of  symmetry,  such  that  Ijy  =  lyz  =  0.  Sideslip 
angle  is  assumed  small  enough  to  use  the  small  angle 
approximation  that  cos  j8  =  1.  Finally,  in  order  to  preserve 
the  linearity  of  the  model,  the  trigonometric  functions  in 
equations  (51)  through  (53)  are  assumed  to  be  constant,  and 
are  evaluated  at  the  initial  conditions.  With  these 
assumptions,  the  equations  become 


&  =  g  + 


gcosBj^cos^^ 

Vcosa^ 


qS 


mVcostti 


6e) 


(54) 
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■J  =  [C*.  ♦  Q.®  *  C’M.Se  *  1 

-Lyy  ^  V 

(55) 

and 

0  =  gcos4>i  +  0(5 

(56) 

The  corresponding  measurement  equations  (7:45)  are 

“■ = 4'‘  - 

(57) 

(58) 

=  0 
in 

(59) 

and 

^  "mg  g  9  g 

(60) 

The  K„  term  in  equation  (57)  corrects  for  boom  bending  and 

upwash.  Since  these  models  were  not  available  for  the 

flight  test  aircraft,  this  term  is  deleted.  The  assumptions 

made  above  also  allow  deletion  of  the  lateral  terms, 

yielding  a  set  of  measurement  equations  of  the  form 

%  =  a  ' 

(61) 
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O'*  =  <3r 


(62) 


e„  =  e 


(63) 


and 


=  (c*,  >  c».«  - 


(64) 


where  the  term  is  expanded  consistent  with  equation  (54)  . 
These  equations  are  used  primarily  to  correct  measured  data 
for  offsets  of  the  respective  sensors  (angle  of  attack  and 
acceleration  packages)  from  the  center  of  gravity  of  the 
aircraft.  In  fact,  the  terms  X„,  and  Z^n  are  those 

offsets  from  the  center  of  gravity  to  either  the  angle  of 
attack  vane  or  the  normal  acceleration  measurement  assembly. 
In  the  Cyber  MMLE3  program,  these  equations  are  incorporated 
into  the  observation  matrices  in  the  model.  However,  to 
ease  processing  required  every  time  the  model  is  accessed, 
the  calculations  are  made  to  correct  the  measured  data  prior 
to  importing  it  into  the  program.  This  simplifies  the 
observation  matrices  in  the  model  significantly. 

If  parameters  for  the  MMLE  model  are  defined  as 


P(l)  -  Cn„ 
P(2)  =  Cm„ 
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P(3) 

~  ^Mq 

P(4) 

P(5) 

~  ^Mie 

the  state  space  model  is  represented  by  the  matrices 


A  = 


■ - — — p(l) 

mVcosa^ 

^P(2) 

^yy 

0 


COs4>i 


B  = 


QS  flrcos6iCos<|)i 

- - p(4)  - 


mVcosa.^ 


qS 


qSc 


p(5) 


yy 


Vcosoi, 


mVcosa^ 


ggC( 

e„ 


-«to 


C  = 


1 

0 

0 


0 

1 

0 


-2^p(l)  0 
mg 


and 


D  = 


0 

0 

0 


0 

0 

0 


2^p(4) 


mg 


mg 


(65) 


(66) 


(67) 


(68) 
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This  is  the  final  form  of  the  model  which  is  coded  for 
using  MULE  on  flight  test  data.  The  model  requires  quite  a 
bit  more  information  to  run  than  the  previous  short  period 
approximation  and  transfer  function  models.  Subsequently,  a 
macro  file  was  created  to  provide  a  user  interface  in  the  PC 
Matlab  environment  to  run  MMLE.  The  resulting  macro  file, 
called  HSINIT.M,  is  shown  in  Figure  7  below. 

The  initialization  macro  file,  when  executed,  prompts 
the  user  to  enter  required  data  interactively.  The  required 
inputs  include  physical  aircraft  parameters,  wind  tunnel 
estimates  of  the  five  stability  derivatives,  and  the 
location  of  the  flight  test  data  file.  In  particular,  the 
required  aircraft  parameters  are 

1.  sref  =  aircraft  reference  area  (S)  in  square  feet. 

2.  cbar  =  mean  aerodynamic  chord  in  feet. 

3.  gw  =  aircraft  gross  weight  in  pounds. 

4.  eg  =  aircraft  center  of  gravity  in  percent  of  mean 
aerodynamic  chord. 

5.  iyy  =  pitch  moment  of  inertia  in  slug-square  feet. 

6.  xab  =  the  X  offset  (positive  forward)  of  the  angle 
of  attack  vane  from  the  center  of  gravity  in  feet. 

7.  xan  =  the  x  offset  (positive  forward)  of  the  normal 
accelerometer  from  the  center  of  gravity  in  feet. 
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clear; 

% 

%  INITIAL  SETUP  MACRO  FOR  RUNNING  HSMMLE.M 
%  MUST  BE  RUN  PRIOR  TO  HSMMLE.M  IN  ORDER  TO 
%  INPUT  REQUIRED  CONSTANTS,  WIND  TUNNEL 
%  ESTIMATES,  AND  LOAD  THE  FLIGHT  TEST  DATA 
% 

% - INPUT  REQUIRED  CONSTANTS 

% 

sref=input( 'REFERENCE  AREA  (S)  IN  SQUARE  FEET?  •); 
cbar=input( 'MEAN  AERODYNAMIC  CHORD  IN  FEET?  • ) ; 
gw=input( 'AIRCRAFT  GROSS  WEIGHT  IN  POUNDS?  '); 
cg=input( 'AIRCRAFT  CG  IN  %MAC?  • ) ; 

iyy=input( 'MOMENT  OF  INERTIA  (lYY)  IN  SLUG-FT"2?  '); 
xab=input( 'X-DIST  IN  FT  FROM  CG  TO  ALPHA  VANE  (+FWD)?  ' ) ; 
xan=input( 'X-DIST  IN  FT  FROM  CG  TO  NZ  ACCEL  {+FWD)?  ' ) ; 
zan=input( 'Z-DIST  IN  FT  FROM  CG  TO  NZ  ACCEL  (+DWN)?  '); 
vtrue=input( 'TRUE  AIRSPEED  IN  FEET  PER  SECOND?  ' ) ; 
qbar=input ( ' DYNAMIC  PRESSURE  IN  LBS  PER  SQUARE  FOOT?  ' ) ; 
% 

% - INPUT  WIND  TUNNEL  ESTIMATES 

% 

pref (l)=input( 'CN_ALPHA  FROM  WIND  TUNNEL  (1/DEG)?  '); 
pref (2)=input( 'CM_ALPHA  FROM  WIND  TUNNEL  (1/DEG)?  ' ) ; 
pref ( 3 ) =input ( ' CM_Q  FROM  WIND  TUNNEL  ( 1/RAD) ?  ' ) ; 
pref (4)=input( 'CN_DE  FROM  WIND  TUNNEL  (1/DEG)?  '); 
pref (5) =input ( ' CM_DE  FROM  WIND  TUNNEL  (1/DEG) ?  ' ) ; 
pref=pref*180/pi; 
pref ( 3 ) =pref ( 3 ) *pi/ 180 ; 

% 

% - LOAD  DATA  FILE 

% 

%  A  FLIGHT  TEST  DATA  FILE  MUST  EXIST  AS  AN  ASCII  FILE 
%  ON  THE  MATLAB  DIRECTORY  (WITH  A  .DAT  EXTENSION) 

% 

data=input ( ' DATA  FILE  (WITHOUT  .DAT  EXTENSION)?  ' , ' s ' ) ; 
ndp=input( 'NUMBER  OF  DATA  POINTS?  '); 
ldc= [ ' load  ' , data , ' . dat ; ' ] ; 
eval^dc)  ; 

% 

dt=l/32; 

t= [ 0 : ndp- 1 ] *dt ; 

% 

save  hsinit 
% 

% - ENDHSINIT.M 


Figure  7.  MMLE  Initialization  Macro  Code 
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8.  zan  -  the  z  offset  (positive  down)  of  the  normal 
accelerometer  from  the  center  of  gravity  in  feet. 

9.  vtrue  =  true  airspeed  in  feet  per  second. 

10.  qbar  =  dynamic  pressure  in  pounds  per  square  foot. 

The  wind  tunnel  values  for  the  stability  derivatives  are 
input  in  common  units  of  per  degree  for  all  but  which 

is  per  radian.  The  macro  then  converts  all  units  to  radians 
to  avoid  dimensional  problems  in  the  model.  The  flight  test 
data  must  exist  in  an  ASCII  file  on  the  PC  Matlab  directory. 
The  macro  is  set  up  to  import  an  ASCII  file  with  a  .dat  file 
extension.  Finally,  the  macro  creates  a  time  vector  with 
the  32  samples  per  second  which  are  available  on  the  flight 
test  aircraft,  as  well  as  the  user  input  number  of  data 
points  in  the  data  sample.  The  initial  data  is  then  saved 
in  a  data  file  to  be  accessed  by  the  MMLE  macro  file. 

The  macro  file  HSP2SS.M  performs  the  parameter-to- 
state-space  conversion.  The  file  is  set  up  to  implement  the 
matrices  shown  in  equations  (65)  through  (68) .  HSP2SS.M  is 
shown  in  Figure  8.  The  macro  calculates  the  bias  terms  in 
these  equations  such  that  the  rates  of  the  state  variables 
are  zeroed,  and  the  load  factor  measurement  is  set  to  its 
initial  value.  The  remainder  of  the  code  is  identical  to 
the  parameter-to-state-space  macros  described  in  the 
simulation  section. 
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f unction [a, phi, gam, c,d,q,xO,dt, rowing, b]=hsp2ss(p) ; 

% - 

%  P2SS  FUNCTION  FOR  HSMMLE.M  -  EXTENSIVELY  COMMENTED 


%  p(l)  =  CN_ALPHA  p(2)  =  CM_ALPHA  p(3)  =  CM_Q 

%  p(4)  =  CN_DE  p(5)  =  CM_DE 

% 

% - PERFORM  INITIAL  CALCULATIONS 

cnstl=(-l*qbar*sref ) / (cos (all) *gw*vtrue/32 . 17) ; 
cnst2=qbar *sref  *cbar / iyy ; 
cnst3=cnst2*cbar/ (2*vtrue) ; 

cnst4=32 . 17*cos (phi) *cos (thl) / (vtrue*cos (all) ) ; 
cnst5=qbar*sref /gw; 

cnst6=-l* (cnstl*p(l) *all+ql+cnstl*p(4) *del+cnst4) ; 
cnst7=-l* (cnst2*p(2) *all+cnst3*p(3) *ql+cnst2*p(5) *del) ; 
cnst8=-l* (ql*cos (phi) ) ; 

cnst9=(-l* (cnst5*p(l) *all+cnst5*p(4) *del) )+nzl; 

% - STABILITY  DERIVATIVES 

a=[cnstl*p(l)  1  0 

cnst2*p(2)  cnst3*p(3)  0 

0  cos (phi)  0]; 

% - CONTROL  DERIVATIVES 

b=[cnstl*p(4)  cnst4+cnst6 

cnst2*p(5)  cnst7 

0  cnstS  ) ; 

% - MEASUREMENT  MATRIX 

C=[  1  0  0 

0  10 
0  0  1 

cnst5*p(l)  0  0  ]; 

% - FEED  THROUGH  MATRIX 

d=[  0  0 

0  0 

0  0 

cnst5*p(4)  cnst9  ]; 

% - STATE  NOISE  COVARIANCE 

q=eye(3)  *le-6;%  SAME  SIZE  AS  A  WITH  Q*Q'  POS.  DEFINITE! 

% - ROWS  IN  Q  IN  WHICH  PARAMETERS  OCCUR 

rowinq=0*p; 

% - INITIAL  STATE  VECTOR 

xO=[all  ql  thl]; 

% - DEFINE  DT  AND  DISCRETIZE 

dt=l/32; 

[phi,gam]=c2d(a,b,dt) ; 

% - end  HSP2SS.M 


Figure  8.  Flight  Test  P2SS  Function  Code 
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All  the  required  up  front  calculations  are  located  in 
the  HSMMLE.M  macro  file  (shown  in  Figure  9) .  The  macro 
first  loads  in  the  file  containing  the  results  of  the 
HSINIT.M  run.  Then,  the  specific  inputs  and  outputs  are 
selected  from  the  flight  test  data  file.  In  this  particular 
case,  there  is  a  correction  of  approximately  four  degrees 
for  a  known  angle  of  attack  calibration  offset.  The  data 
are  also  converted  to  radians  for  unit  consistency  in  the 
program.  Angle  of  attack  and  load  factor  measurements  are 
corrected  for  the  distance  of  the  sensors  from  the  center  of 
gravity,  and  initial  values  for  the  states,  measurements, 
and  control  inputs  are  calculated  from  the  flight  test  data 
file.  The  parameter  identification  options  are  then  set. 

The  only  difference  in  these  setting  from  the  simulation  is 
that  the  initial  parameter  guess  is  set  equal  to  the  wind 
tunnel  estimates.  The  innovations  covariance  matrix  is  set 
to  the  default  for  MMLE3,  which  has  been  developed  over  the 
years  as  a  good  rule  of  thumb  for  flight  test  data  (5:12). 
The  last  section  of  the  macro  file  performs  final 
calculations  to  convert  the  stability  derivative  estimates 
back  into  common  units,  calculate  equivalent  second-order 
system  parameters,  and  display  the  outputs.  The  second- 
order  system  damping  ratio  and  natural  frequency  are 
calculated  using  the  final  system  A  matrix  and  the  PC  Matlab 
DAMP  function.  The  DAMP  function  returns  these  parameters 
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format  compact, clc 
clear; 

load  hsinit; 

! erase  hsmmle . log ; 
diary  hsmmle. log 

I _ _ _ _ _ I . 

dispCHAVE  STABILITY  PARAMETER  IDENTIFICATION  MACRO  FILE') 
dispC  T-38A  LONGITUDINAL  STAB.  AND  CONTROL  DERIVATIVES') 
dispC  ONE  INPUT,  FOUR  OUTPUTS') 

global  sref  cbar  gw  iyy  vtrue  qbar  all  thl  phi  ql  nzl  del; 
% 

% - INPUT  FLIGHT  TEST  DATA 

uydata=zeros (ndp, 6) ; 

coll=[ 'uydata (;,!)=' ,data, * ( : ,7) ; ' ] ; 
eval (coll) ; 

uydata ( : , 1) =uydata( : , 1) *pi/180; 
uydata ( : , 2 ) =ones ( ndp , 1 ) ; 

col3=( ' uydata ( : , 3)=' ,data, ' ( : , 9) +(4 . 07*ones(ndp, 1) ) ; ' ] ; 
eval(col3) ; 

uydata (: ,3)=uydata(: ,3)*pi/180; 
col4= [ ' uydata ( ; , 4 ) = ' , data , ' ( : , 8 ) ; ' ] ; 
eval (col4) ; 

uydata ( ; ,4)=uydata( : ,4) *pi/180; 
col5=( 'uydata( ; , 5)  =  ' ,data, ' ( : , 10)  ; '  ] ; 
eval (col5) ; 

uydata ( : ,5)=uydata( : ,5) *pi/180; 
col6=[ ' uydata ( ; , 6) =' ,data, ' ( : , 12) ; ' ] ; 
eval (cole ) ; 

% 

%  CORRECT  ALPHA  AND  NZ  FOR  DISTANCE  OF  SENSORS  FROM  CG 
uydata ( : , 3 ) =uydata ( : , 3 ) + (xab*uydata ( ; , 4 ) /vtrue) ; 
for  i=l:ndp, 

q2t(i)=uydata(i,4) *uydata(i,4) ; 
q2=q2t ' ; 

end 

uydata ( : , 6)=uydata( : , 6) - (zan*q2/32 . 17) ; 

qdot=[diff (uydata( : ,4) ) *32.0;0] ; 

uydata ( : , 6 ) =uydata ( : , 6 ) - ( xan*qdot / 32.17) ; 

%  INITIAL  CONDITIONS  FOR  ALPHA,  THETA,  PHI,  Q,  NZ,  AND  DE 
a 1 l=uydata (1,3) ; 
thl=uydata (1,5) ; 

phlic=[ 'phl=' ,data, ' (1, 11) *pi/180; ' ] ; 

eval (phlic) ; 

ql=uydata (1,4) ; 

nz l=uydata (1,6) ; 

del=uydata (1,1) ; 


Figure  9.  Flight  Test  MMLE  Macro  Code  (Sheet  1  of  2) 
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% - SUFFICIENT  INPUTS  TO  MMLE  FOLLOW 

pidq=[l] ;%  IDENTIFY  WHICH  PARAMETERS  ARE  TO  BE  IDENTIFIED 
pidm=[l];%  IN  THE  QUADRATIC,  MARQUARDT,  AND  FINAL  STAGES. 
pidf=[l:5] ;%pidq,  pidm,  pidf  MUST  BE  VALID  IF  NOT  USED 
pO=pref;%  INITIAL  PARAMETER  GUESS  =  WIND  TUNNEL  ESTIMATES 
opt=[0  0  10  0  .02  .001  .001  1];%  MAXIMUM  ITERATIONS  AND 

%  CONVERGENCE  CRITERIA 

radc=180/pi ; 

ggiv=[10*(radc"2)  60*(radc"2)  30*(radc"2)  200]; 
ggO=inv(diag(ggiv) ) ;%  INNOVATIONS  COV.  (MMLE3  DEFAULT) 
p2snain=‘hsp2ss' ;  %  FILENAME  OF  P2SS  FUNCTION 

pert=le-4;  %  PERTURBATION  USED  FOR  GRADIENT  CALCULATION 
1 inesearch=l ; 
ininle 

% - PERFORM  FINAL  CALCULATIONS 

cna=pfin(l) /radc; 

cina=pf  in  ( 2 )  /radc ; 

cinq=pfin(3)  ; 

cnde=pf in ( 4 ) / radc ; 

cinde=pf  in  ( 5 )  /radc ; 

constf=sqrt (qbar*sref *cbar/iyy) ; 

[  wn ,  zeta  ]  =dainp  ( a )  ; 

ktheta=cmde*radc*qbar*sref *cbar / iyy ; 
invth2=cna*radc*32 . 17*qbar*sref / (vtrue*gw) ; 
disp('  ') 

dispC  DATA  FILE  =  ' )  ,disp{data) 
disp('  ') 

cond=[cg  all*180/pi  qbar] ; 

dispC  CG  TRIM  AOA  QBAR') 

disp('  ') 

disp(cond) 

disp{'  ') 

disp('  ') 

deriv=[cna  cma  cmq  cnde  cmde] ; 

dispC  CNA  CMA  CMQ  CNDE  CMDE') 

dispC  ') 

disp(deriv) 

dispC  ') 

dispC  ') 

systein=[zeta(3)  wn(3)  ktheta  invth2); 

dispC  ZETA  WN  KTHETA  1/T_THETA2 ' ) 

dispC  ') 

disp( system) 

diary  off 

! print  hsmmle.log; 

% - ENDHSMMLE.M 


Figure  9.  Flight  Test  MMLE  Macro  Code  (Sheet  2  of  2) 
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for  the  input  A  matrix.  The  system  parameters  and  1/T^ 
are  calculated  as 


and 


(69) 


(70) 


The  end  result  of  the  HSINIT.M,  HSP2SS.M,  and  HSMMLE.M  macro 
files  is  a  user-friendly,  interactive  capability  to  perform 
parameter  estimation  from  flight  test  data  using  a  personal 
computer.  The  ultimate  test  of  the  system,  however,  is  to 
accomplish  a  flight  test  program  to  verify  the  results.  The 
following  section  will  describe  the  details  of  the  flight 
test  planning,  execution,  and  results. 
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V.  Flight  Test 


Introduction 

This  section  presents  the  results  of  a  limited 
development  and  evaluation  of  a  parameter  estimation 
capability  for  the  USAF  Test  Pilot  School  (USAFTPS) ,  using 
the  final  model  presented  in  the  last  section.  To  assess 
this  capability,  flight  test  data  were  gathered  in  a  T-38A 
aircraft. 

Several  past  attempts  have  been  made  to  demonstrate 
parameter  estimation  capability  at  the  USAFTPS.  The  most 
recent  was  the  T-38  Parameter  Estimation  Study  (HAVE  PEST) , 
which  was  conducted  from  September  through  December  of  1989 
These  past  projects  suffered  from  problems  with  both  the 
data  acquisition  system  and  computer  software  availability. 
Appropriate  steps  were  taken  to  resolve  these  problems  for 
this  test. 

The  test  team  flew  9  test  sorties  and  3  support 
sorties,  all  in  the  T-38A,  for  a  total  of  12.9  flight  hours 
The  test  missions  were  flown  from  17  September  through  19 
October  1990  at  the  AFFTC,  Edwards  AFB,  California. 

The  overall  test  objective  was  to  evaluate  the 
parameter  estimation  capability  for  the  USAFTPS, 
concentrating  only  on  the  longitudinal  modes  of  a  T-38A 
aircraft. 
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The  specific  objectives  were  to: 

1.  Evaluate  the  personal  computer  based  approach  to 
aircraft  parameter  estimation  using  the  MMLE  routine  in  PC 
Matlab.  Specifically,  to  demonstrate  the  capability  to 
process  data  from  the  USAFTPS  data  reduction  facility 
through  these  programs  to  accomplish  the  parameter 
estimation  task. 

2 .  Compare  results  obtained  using  this  method 
available  at  USAFTPS  with  results  of  the  MMLE3  program  run 
on  the  Ridley  Mission  Control  Center  (RMCC)  Cyber  computer. 

3.  Compare  equivalent  system  parameters  obtained  from 
the  MMLE  results  to  those  obtained  from  frequency  response 
analysis  (FRA)  methods  available  on  the  RMCC  Cyber  computer, 
as  well  as  from  classical  modal  analysis  flight  test 
techniques. 

Test  Item  Description 

The  personal  computer  hardware  used  for  this  test  was  a 
Dell  386/25  with  a  387  Math  Coprocessor.  The  primary 
software  used  was  the  PC  Matlab  program. 

The  PC  Matlab  program  was  a  matrix  operation  program 
which  could  be  used  to  solve  practical  engineering  and 
mathematical  problems.  PC  Matlab  was  designed  for  specific 
applications  such  as  control  theory  and  signal  processing. 
(8:3).  Two  additional  "toolboxes”  were  used  for  this 
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project:  the  Control  System  Toolbox  (6)  and  the  State  Space 

Identification  Tool  (10)  . 

Flight  test  data  from  the  T-38A  were  used  as  a  tool  to 
evaluate  the  parameter  estimation  capabilities  of  the 
USAFTPS.  The  test  aircraft,  USAF  serial  number  (S/N)  63- 
8135,  was  a  production  aircraft  modified  with  the  flight 
test  instrumentation  described  below.  The  T-38A,  built  by 
Northrop  Corporation,  was  a  supersonic,  two-place,  tandem 
trainer  powered  by  two  afterburning  J85-GE-5  series  turbojet 
engines.  The  flight  controls  were  powered  by  a 
conventional,  irreversible  hydraulic  system.  The  aircraft 
was  flown  with  no  external  stores.  More  information  about 
the  production  model  T-38A  can  be  found  in  the  T-38  Flight 
Manual  (17) . 

The  test  aircraft  was  modified  as  follows: 

1.  Nose  boom  with  a  Yaw,  Angle-of -attack,  Pitot- 
static  System  (YAPS)  head. 

2.  Metraplex  airborne  test  instrumentation  system 
with  32  samples  per  second  data  collection 
rate. 

3.  Sensitive  airspeed,  Mach,  and  G  indicators. 

4 .  C-band  beacon . 

Instrumentation 

Table  A1  in  Appendix  A  shows  the  instrumentation 
parameters  and  available  resolutions  for  this  test.  The 
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Metraplex  system  was  essential  to  the  success  of  the  test 
program.  Its  32  sample  per  second  sampling  rate  and  greater 
parameter  resolution  were  a  must  for  the  parameter 
estimation  process.  In  general,  instrumentation  issues  were 
critical  to  achieving  test  objectives.  Several  key  lessons 
were  learned  regarding  instrumentation,  and  these  are 
summarized  below. 

Calibrations.  A  key  to  the  success  of  the  program  was 
ensuring  that  the  instrumentation  was  properly  calibrated. 

To  accomplish  this,  the  test  team  developed  specific 
calibration  requirements  (including  calibration  ranges  and 
increments) .  To  ensure  the  quality  of  the  calibrations, 
test  team  members  participated  actively  in  the  entire 
calibration  process.  Critical  parameters  (angle  of  attack, 
pitch  rate,  normal  acceleration,  and  elevator  deflection) 
were  calibrated  three  times  to  minimize  errors.  In  some 
cases,  the  calibration  curves  were  nonlinear  at  the  extremes 
of  the  calibration  ranges.  The  USAFTPS  data  reduction 
system  could  only  implement  first  order  (i.e.  linear)  curve 
fits.  In  those  instances  where  nonlinear  calibrations 
existed,  the  ranges  of  the  calibrations  were  reduced  to 
accurately  represent  the  regions  of  primary  interest. 

Parameter  Resolution.  Resolution  of  key  parameters  was 
also  critical  to  the  success  of  the  parameter  estimation 
effort.  Based  on  the  data  acquisition  system  configuration, 
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calibration  ranges  were  narrowed  to  achieve  desired 
resolutions.  Even  so,  the  pitch  angle  resolution  turned  out 
to  be  inadequate.  Better  resolution  could  have  been 
achieved  by  narrowing  the  calibration  range  even  further. 

Parameter  Filtering.  Data  from  four  otherwise 
effective  sorties  were  corrupted  by  a  filter  mismatch  among 
the  key  longitudinal  parameters.  The  data  acquisition 
system  was  set  up  to  filter  some  parameters  with  10  Hertz 
filters,  some  with  20  Hertz  filters,  and  some  not  at  all. 
This  created  a  problem  with  frequency  response  analysis 
techniques,  and  to  a  limited  degree  with  the  MMLE 
techniques.  Filtering  of  key  parameters  should  be  checked 
for  consistency  to  avoid  inducing  unwanted  phase  lag 
differences  between  parameters. 

Other  Issues.  Several  inaccuracies  were  introduced 
which  caused  instrumentation  measurements  to  differ  from 
true  aircraft  response.  The  positions  of  accelerometer 
packages  and  nose  boom  measurement  vanes  relative  to 
aircraft  center  of  gravity  directly  affected  acceleration 
and  angle  measurements  from  those  sources.  In  fact, 
vertical  displacements  could  not  be  determined  due  to 
uncertainty  of  the  vertical  datum  for  the  aircraft. 

Accurate  measurements  of  these  relative  positions  are  a  must 
for  getting  accurate  parameter  estimation.  In  addition, 
corrections  could  have  been  made  for  boom  bending  and  upwash 
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in  angle  of  attack  measurements.  No  such  models  were 
readily  available,  and  these  corrections  were  not  made  for 
this  test.  More  accurate  parameter  estimates  could  be 
obtained  by  including  such  a  model. 

Test  and  Evaluation 

Test  Procedures.  Data  were  collected  for 
classical  analysis  using  pitch  doublets  at  each  test  point. 
Figure  A1  in  Appendix  A  shows  the  four  test  points  and  their 
positions  in  the  T-38A  operational  flight  envelope.  Table 
A2  describes  each  test  point  in  detail.  Doublets  were 
performed  at  approximately  the  short  period  natural 
frequency  of  the  aircraft.  These  maneuvers  and  the 
corresponding  hand  reduction  of  the  data  are  described  in 
the  Flying  Qualities  Phase  Textbook  (1:8.47-8.53). 

Flight  test  data  were  gathered  for  MMLE  techniques 
using  pitch  doublets  at  each  test  point.  Doublets  were 
repeated  at  elevated  angles  of  attack  (in  increments  to  give 
at  least  two  test  points  above  trim  without  encountering 
buffet) .  Load  factor  was  increased  to  provide  the  increase 
in  angle  of  attack,  with  altitude  loss  accepted  to  maintain 
Mach  number.  Doublets  were  performed  with  a  period  of  one 
second.  Amplitude  was  adjusted  to  get  a  pitch  rate  of 
approximately  ±10  degrees  per  second.  Angle  of  attack 
perturbations  were  kept  within  ±3  degrees.  These  maneuvers 
are  described  in  Flight  Investigation  of  Various  Control 
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Inputs  Intended  for  Parameter  Estimation  (15) .  Each  pilot 
practiced  the  MMLE  inputs  in  a  ground  simulator  prior  to  the 
start  of  flight  test.  Data  were  processed  through  the  PC 
Mat lab  and  Cyber  MMLE  programs,  which  were  discussed  in 
detail  in  the  previous  sections.  The  programs  were  used  to 
generate  stability  derivatives,  as  well  as  estimated  second- 
order  system  parameters  based  on  these  derivatives.  The 
required  aircraft  gross  weight  and  moment  of  inertia  data 
were  obtained  from  Calculated  Dead  Weight  Distribution  and 
Moment  of  Inertia  Report.  T-38A  (13) .  The  correlation  of 
this  data  to  the  actual  flight  test  aircraft  was  in 
question,  but  it  was  the  best  data  available.  Future 
attempts  to  perform  parameter  estimation  on  the  T-38A  should 
calculate  these  data  more  precisely  for  the  test  aircraft. 
Wind  tunnel  estimates  were  used  as  starting  points  for  the 
derivative  estimation  processes,  and  were  obtained  from  T-38 
Static  Stability  and  Basic  Aerodynamic  Data  (2) . 

Pitch  frequency  sweeps  were  used  to  generate  data  for 
frequency  response  analysis.  Frequency  sweeps  were 
performed  from  approximately  0.2  to  20  radians  per  second, 
with  the  entire  sweep  taking  at  least  32  seconds.  Each 
pilot  practiced  the  frequency  sweep  inputs  in  a  ground 
simulator  prior  to  the  start  of  flight  test.  Data  were 
processed  through  the  Cyber  FRA  program.  The  FRA  program 
converts  time  history  data  into  the  frequency  domain  using  a 
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Fourier  transform,  and  produces  output  frequency,  magnitude, 
and  phase  angle  information  (16:826).  This  information  was 
matched  using  an  equivalent  system  curve  fit,  in  order  to 
obtain  equivalent  second-order  system  parameters.  The 
results  were  used  to  compare  against  the  same  parameters 
calculated  from  the  MMLE  stability  derivative  estimates. 

This  provided  a  means  of  checking  the  MMLE  results  using  a 
completely  different  (frequency  domain)  technique. 

The  frequency  sweeps  were  difficult  to  perform  well. 

The  technique  required  a  very  low  frequency  input  initially 
(about  0.2  radians  per  second).  Pilots  found  it  challenging 
to  make  such  a  slow  input  and  keep  the  amplitude  small 
enough  to  avoid  large  excursions  in  altitude  and  airspeed. 

A  smooth  increase  in  frequency  was  then  required,  with  input 
amplitude  adjusted  to  keep  essentially  a  constant  output 
amplitude.  At  the  higher  frequencies,  pilots  had  to  use 
some  outside  reference  to  decouple  themselves  from  the 
response  of  the  aircraft.  Practice  helped  the  quality  of 
the  inputs  significantly,  and  pilots  generally  felt  that  the 
time  spent  in  the  ground  simulator  prior  to  flight  was 
invaluable.  Ground  simulator  training  should  be 
accomplished  before  collecting  inflight  data  using  parameter 
estimation  flight  test  techniques. 

Test  Results.  The  data  from  the  MMLE  doublets  were 
successfully  processed  on  the  PC.  Flight  test  data  were 
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reduced  to  engineering  units  through  the  Test  Pilot  School 
PDP  11/84  computer,  as  detailed  in  Data  Reduction  on  the 
Test  Pilot  School  Data  Reduction  Computer  System  (3) .  The 
existing  plot  building  routine  was  used  to  create  data  files 
with  the  specific  time  slices  for  MMLE  processing.  These 
files  could  easily  be  downloaded  to  the  PC  using  the  file 
transfer  protocol  (FTP)  through  the  existing  USAFTPS 
personal  computer  network.  The  PC  Matlab  program  was 
capable  of  reading  in  ASCII  files  such  as  those  created  by 
the  plot  building  program,  provided  the  header  lines  had 
been  deleted. 

Although  the  processing  procedures  for  MMLE  on  the  PC 
were  quite  lengthy,  they  were  easily  executed  using  the 
macro  files  discussed  in  previous  sections.  The  MMLE  macro 
files  were  set  up  to  be  as  generic  as  possible  and  allow  the 
user  to  input  all  required  program  inputs  interactively  at 
the  terminal. 

Figure  10  below  shows  an  example  of  the  time  history 
matching  by  PC  MMLE  for  test  point  1(a).  Figures  B1  through 
Bll  in  Appendix  B  show  the  remaining  time  history  matching 
results  from  the  PC  MMLE.  In  general,  the  lack  of  boom 
bending  and  upwash  corrections  to  angle  of  attack  can  be 
seen  in  the  lower  quality  of  the  match  to  the  flight  test 
angle  of  attack  data.  In  two  particular  cases.  Figures  B4 
and  B8,  the  flight  test  data  for  angle  of  attack  were 
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Pitch  Rate,  q  (radians  /  second)  Angle  of  Attack,  a  (radians) 


T-38A  USAF  S/N  63-8135 
CRUISE  CONFIGURATION 
NO  EXTERNAL  STORES 


GROSS  WEIGHT  -  10,258  LBS 

CG:  16.6  %  MAC 

YAW  AUGMENTATION:  OFF 


Figure  10.  Time  History  Match,  Test  Point  1(a) 

(Sheet  1  of  2) 
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Ljoad  Factor,  Nz  (g)  Pitch  Angle,  6  (radians) 


T-38A  USAF  S/N  63-8135 
CRUISE  CONFIGURATION 
NO  EXTERNAL  STORES 


GROSS  WEIGHT  -  10,258  LBS 

CG:  16.6  %  MAC 

YAW  AUGMENTATION:  OFF 


Figure  10.  Time  History  Match,  Test  Point  1(a) 

(Sheet  2  of  2) 
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erroneous.  However,  the  program  was  still  able  to  complete 
the  parameter  identification  process  with  reasonable 
results. 

Figures  11  through  15  on  the  following  pages  show  the 
comparison  of  the  results  from  the  PC  and  the  Cyber  for  MMLE 
estimation  of  longitudinal  stability  derivatives.  The 
derivatives  were  estimated  based  on  the  Test  Pilot  School 
sign  convention,  as  shown  in  Figure  A2  in  Appendix  A. 
Different  symbols  are  presented  for  the  PC  and  Cyber  results 
at  the  two  dynamic  pressures  tested.  The  vertical  bars 
indicate  the  modified  Cram6r-Rao  bounds  on  the  Cyber 
estimates.  The  Cramer-Rao  bounds  are  multiplied  by  a  factor 
of  ten,  which  is  used  to  indicate  the  uncertainty  level  in 
the  parameter  estimates.  Recall  that  the  Cram6r-Rao  bounds 
are  the  lower  limit  on  the  variance  of  the  estimates,  and 
are  frequently  multiplied  by  a  scale  factor  to  arrive  at 
uncertainty  levels  (4:11).  These  bounds  are  used  in  lieu  of 
confidence  intervals,  and  indicate  the  confidence  in  the 
Cyber  estimates  for  each  data  run.  A  detailed  description 
of  Cram6r-Rao  bounds  can  be  found  in  the  User's  Manual  for 
MMLE3  (7) .  The  results  from  the  PC  correlated  well  with 
those  from  the  Cyber.  In  all  cases,  the  PC  estimates  were 
within  the  uncertainty  bounds.  In  fact,  the  vast  majority 
of  the  PC  runs  fell  well  within  the  bounds  of  the  Cyber 
estimates.  Overall,  the  quality  of  the  PC  MMLE  results 
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T-38A  USAF  S/N  63-81 35  AVERAGE  WEIGHT:  1 0,571  LBS 

CRUISE  CONRGURATION  AVERAGE  CG:  17.1  %  MAC 
NO  EXTERNAL  STORES  YAW  AUGMENTATION:  OFF 


Angle  of  Attack  (degrees) 


T-38A  USAF  S/N  63-8135  AVERAGE  WEIGHT:  10,571  LBS 

CRUISE  CONFIGURATION  AVERAGE  CG:  17.1  %  MAC 
NO  EXTERNAL  STORES  YAW  AUGMENTATION:  OFF 


71 


Angle  of  Attack  (degrees) 


T-38A  USAF  S/N  63-81 35  AVERAGE  WEIGHT:  1 0,571  LBS 

CRUISE  CONRGURATION  AVERAGE  CG:  17.1  %  MAC 
NO  EXTERNAL  STORES  YAW  AUGMENTATION:  OFF 
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Angle  of  Attack  (degrees) 


T-38A  USAF  S/N  63-8135  AVERAGE  WEIGHT:  10,571  LBS 

CRUISE  CONFIGURATION  AVERAGE  CQ:  17.1  %  MAC 
NO  EXTERNAL  STORES  YAW  AUGMENTATION:  OFF 
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Figure  14.  MMLE  Simmary 


T-38A  USAF  S/N  63-8135  AVERAGE  WEIGHT:  10,571  LBS 

CRUISE  CONFIGURATION  AVERAGE  CG:  17.1  %  MAC 
NO  EXTERNAL  STORES  YAW  AUGMENTATION:  OFF 
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obtained  was  satisfactory  for  integration  into  the  Test 
Pilot  School  curriculum.  However,  the  scope  of  the 
evaluation  was  limited  to  longitudinal  aircraft  responses. 
The  capability  of  the  Test  Pilot  School  to  identify  lateral- 
directional  responses  should  be  evaluated  prior  to 
incorporating  MMLE  into  the  flying  qualities  course. 

Figure  16  below  shows  an  example  of  the  Cyber  FRA 
results  from  test  point  1(a).  Figures  B12  through  B14  in 
Appendix  B  present  the  remaining  Cyber  results.  Data  from 
multiple  time  histories  were  ensemble  averaged  to  improve 
the  quality  of  the  estimation.  This  project  achieved 
excellent  results  from  FRA  using  frequency  sweep  inputs. 

The  use  of  ensemble  averaging  was  a  key  factor  in  generating 
such  high  quality  results.  Ensemble  averaging  should  be 
considered  a  must  for  future  attempts  at  frequency  response 
analysis.  The  data  are  presented  in  Bode  plot  form,  with 
magnitude  and  phase  angle  of  the  pitch  rate  to  elevator 
deflection  transfer  function  plotted  versus  frequency.  The 
Cyber  plots  show  estimates  with  95  percent  confidence 
intervals.  In  general,  the  Cyber  results  were  excellent. 

The  confidence  intervals  were  tight  out  to  approximately  15 
radians  per  second.  A  lower-order  equivalent  system  (LOES) 
curve  fit  was  used  to  match  the  frequency  characteristics 
from  Figures  16  and  B12  through  B14 .  The  matching  results 
are  plotted  in  Figures  17  through  20  on  the  following  pages. 
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T-38AUSAFS/N  63-81 35  AVERAGE  WEIGHT:  10,105  LBS 

CRUISE  CONFIGURATION  AVERAGE  CQ:  16.7  %  MAC 
NO  EXTERNAL  STORES  YAW  AUGMENTATION:  OFF 


LEGEND 
o  Cyber  Estimates 

20-  -  9516  Confidence  Limits 


-O 


■“oil  '  r—^—r-T-r^  ,  '  '  '  '  "I'o - ' - '  '  '  '  '  ’  'll 

Frequency  (radians  /  second) 


Frequency  (radians  /  second) 


Figure  16.  FRA  Results,  Test  Point  1(a) 
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T-38A  USAF  S/N  63-8135 
CRUISE  CONFIGURATION 
NO  EXTERNAL  STORES 


AVERAGE  WEIGHT:  10,105  LBS 
AVERAGE  CG:  16.7  %  MAC 
YAW  AUGMENTATION:  OFF 


L£GEND 
Right  Test  Data 
Lower  Order  System  Rt 


Frequency  (radians  /  second) 


Frequency  (radians  /  second) 


Figure  17.  LOSS  Natch,  Test  Point  1(a) 
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Figure  18.  LOES  Match,  Test  Point  2(a) 
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Figure  19.  LOES  Natch,  Test  Point  3(a) 
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Figur*  20.  LOSS  Match,  Test  Point  4(a) 
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Finally,  a  comparison  of  equivalent  second-order  system 
parameters  is  presented  in  Tables  5  through  8  below  for  each 


test  point.  Classical  analysis  techniques  were  used  to 
calculate  damping  ratio  (f^)  and  natural  frequency  (o)^)  from 
the  classical  pitch  doublet  inputs  (1:8.48-8.53).  The 
stability  derivatives  estimated  by  the  PC  and  Cyber  MMLE 
runs  were  used  to  approximate  all  four  system  parameters 
shown,  as  described  in  the  last  section.  The  LOES  match 
discussed  above  was  used  to  determine  the  system  parameters 
from  Cyber  FRA. 

In  most  cases,  the  PC  results  for  each  test  point  were 
less  than  five  percent  different  from  the  Cyber.  In  all 
cases,  the  PC  results  were  within  10  percent.  Overall,  the 
system  parameters  correlated  extremely  well  among  the 
different  data  collection  methods  and  data  reduction 
techniques. 


Table  5.  System  Parameter  Summary,  Test  Point  1(a) 


MMLE 

MMLE 

FRA/ LOES 

PARAMETER 

CLASSICAL 

(PC) 

(CYBER) 

(CYBER) 

fn. 

0.31 

0.3584 

0.3475 

0.3180 

a)„  (rad/s) 

2.54 

2.8599 

2.8241 

2.9680 

K, 

— 

16.4185 

15.9954 

15.3100 

1/T« 

— 

1.1080 

1.1852 

1.2470 

Table  6.  system  Parameter  Summary,  Test  Point  2(a) 


MMLE 

MMLE 

FRA/LOES 

PARAMETER 

CLASSICAL 

(PC) 

(CYBER) 

(CYBER) 

0.32 

0.3412 

0.3275 

0.3103 

w«,  (rad/s) 

3.26 

4.0590 

4.0193 

3.9680 

K, 

— 

29.4538 

28.4696 

25.6100 

1/T« 

— 

1.4363 

1.5946 

1.6680 

Table  7.  System  Parameter  Summary,  Test  Point  3(a) 


MMLE 

MMLE 

FRA/LOES 

PARAMETER 

CLASSICAL 

(PC) 

(CYBER) 

(CYBER) 

U 

0.31 

0.2872 

0.2703 

0.3196 

(rad/s) 

2.64 

2.8034 

2.7365 

2.7820 

K, 

— 

18.4411 

17.8397 

18.4400 

1/T« 

— 

0.8924 

0.8505 

0.7041 

Table  8.  System  Parameter  Summary,  Test  Point  4(a) 


MMLE 

MMLE 

FRA/LOES 

PARAMETER 

CLASSICAL 

(PC) 

(CYBER) 

(CYBER) 

0.29 

0.3079 

0.3244 

0.2880 

aj„  (rad/s) 

3.28 

3.8058 

3.8660 

3.9520 

K, 

— 

32.0179 

33.4768 

30.6200 

1/T« 

— 

1.2477 

1.1767 

1.3750 

VI.  Conclusions  and  Recommendations 


A  system  of  macro  files  has  been  developed  for 
performing  Modified  Maximum  Likelihood  Estimation  of 
aircraft  stability  derivatives  using  a  personal  computer. 

The  process  is  completed  in  the  environment  of  the  PC  Matlab 
personal  computer  program,  with  integration  of  flight  test 
data  files  from  the  Test  Pilot  School  data  reduction 
facility. 

Results  from  this  personal  computer  based  approach  have 
been  compared  with  results  from  the  mainframe  computer 
program  (MMLE3)  which  is  considered  the  standard  for 
parameter  estimation  at  the  Air  Force  Flight  Test  Center. 

The  personal  computer  results  were  generally  well  within  the 
uncertainty  levels  of  the  mainframe  estimates.  A  comparison 
of  second-order  system  parameters  obtained  from  classical 
analysis,  stability  derivatives  extracted  using  both  MMLE 
programs,  and  frequency  response  analysis  techniques  showed 
that  the  system  parameters  correlated  extremely  well  among 
the  different  data  collection  methods  and  data  reduction 
techniques.  Overall,  the  accuracy  of  the  personal  computer 
method  was  satisfactory,  verifying  the  feasibility  of  the 
personal  computer  based  approach  to  aircraft  parameter 
identification . 
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The  quality  of  flight  test  inputs  was  critical  to 
successful  parameter  estimation.  Frequency  sweeps  were 
particularly  challenging  to  perform.  Practice  helped  the 
quality  of  the  inputs  significantly,  and  pilots  felt  that 
the  time  spent  in  the  ground  simulator  prior  to  flight  was 
invaluable.  Ground  simulator  training  should  be 
accomplished  before  collecting  inflight  data  using  parameter 
estimation  flight  test  techniques. 

The  best  data  available  for  moments  of  inertia  were  in 
a  Northrop  Corporation  report  from  1958.  Future  attempts  to 
perform  parameter  estimation  on  the  T-38A  should  calculate 
data  which  more  closely  represents  the  flight  test  aircraft. 
This  would  ensure  estimates  are  as  close  as  possible  to 
those  of  the  actual  aircraft. 

Other  possibilities  for  future  research  projects 
include  investigation  of  including  boom  bending  and  upwash 
corrections  to  the  angle  of  attack  measurements.  This  would 
improve  the  time  history  match  of  the  angle  of  attack  trace. 
Another  possible  area  for  improvement  is  correction  for  non¬ 
zero  lateral  and  directional  parameters.  The  current  system 
assumes  that  for  a  longitudinal  pulse  there  are  no  lateral 
or  directional  responses.  It  would  be  an  improvement  to  the 
flexibility  of  the  model  to  incorporate  these  corrections 
into  the  model.  Finally,  PC  MMLE  capability  for  the 
lateral-directional  axes  could  be  investigated. 
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Appendix  A;  Flight  Test  Information 


Table  Al.  Test  Parameter  Ranges  and  Resolutions 


RESOLUTION 


PARAMETER 


TIME 


ELEVATOR  POSITION 


LEFT  AILERON  POSITION 


RIGHT  AILERON  POSITION 


RUDDER  POSITION 


LONG.  STICK  POSITION 


LATERAL  STICK  POSITION 


RUDDER  PEDAL  POSITION 


LONG.  STICK  FORCE 


LATERAL  STICK  FORCE 


RUDDER  PEDAL  FORCE 


ANGLE  OF  ATTACK 


ANGLE  OF  SIDESLIP 


PITCH  ANGLE 


BANK  ANGLE 


PITCH  RATE 


ROLL  RATE 


YAW  RATE 


NORMAL  ACCELERATION 


LONG.  ACCELERATION 


LATERAL  ACCELERATION 


RANGE 


-8  to  +16  deg 


-35  to  +25  deg 


-25  to  +35  deg 


±30  deg 


-5  to  +7  in 


±8  in 


±5  in 


±70  lb 


±35  lb 


±150  lb 


-10  to  +40  deg 


±10  deg 


±60  deg 


±90  deg 


±50  deg/s 


±100  deg/s 


±50  deg/s 


-2  to  +5  g 


-1  to  +1  g 


-1  to  +1  g 


millisecond 

0.03 

deg 

0.06 

deg 

0.06 

deg 

0.06 

deg 

0.02 

in 

0.02 

in 

0.01 

in 

0.20 

lb 

0.10 

lb 

0.20 

lb 

0.05 

deg 

0.02 

deg 

0.20 

deg 

0.40 

deg 

0.10  deg/s 

0.20  deg/s 

0.10  deg/s 

0.01 

g 

0.01 

g 

0.01 

g 
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Figure  Al.  Flight  Test  Envelope 


Table  A2.  Test  Point  Description 


TEST  POINT 

PRESSURE 

ALTITUDE 

MACH 

NUMBER 

DYNAMIC 

PRESSURE 

10,000  ft 

0.50 

250  psf 

2  (a),  (b) ,  and  (c) 

10,000  ft 

0.65 

430  psf 

3  (a),  (b) ,  and  (c) 

23,000  ft 

0.65 

250  psf 

4  (a)  ,  (b)  ,  and  (c) 

23,000  ft 

0.85 

430  psf 

Note:  At  each  test  point,  the  (a) ,  (b) ,  and  (c) 

correspond  to  angles  of  attack  of  Ig  trim,  trim  +  2 
degrees,  and  trim  +  4  degrees,  respectively. 
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Figure  A2.  Test  Pilot  School  Sign  Convention 
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Figure  Bl.  Time  History  Match,  Test  Point  1(b) 

(Sheet  1  of  2) 

8B 


Figure  Bl.  Time  History  Match,  Test  Point  1(b) 

(Sheet  2  of  2) 
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Pitch  Rate,  q  (radians  /  second)  Angle  of  Attack,  a  (radians) 


T-38A  USAF  SlU  63-81 35  GROSS  WEIGHT  -  9,968  LBS 

CRUISE  CONFIGURATION  CG;  1 6.8  %  MAC 
NO  EXTERNAL  STORES  YAW  AUGMENTATION:  OFF 


Time  (seconds) 


Time  (seconds) 


Figure  B2.  Time  History  Match,  Test  Point  1(c) 

(Sheet  1  of  2) 
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Load  Factor,  Nz  (g)  Pitch  Angle,  6  (radians) 


Figure  B2.  Time  History  Match^  Test  Point  1(c) 

(Sheet  2  of  2) 
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Pitch  Rate,  q  (radians  /  second)  Angle  of  Attack,  a  (radians) 


T-38AUSAFS/N  63-81 35  GROSS  WEIGHT- 10,378  LBS 

CRUISE  CONFIGURATION  CG:  16.5  %  MAC 
NO  EXTERNAL  STORES  YAW  AUGMENTATION:  OFF 


Time  (seconds) 


Figure  B3.  Time  History  Natch,  Test  Point  2(a) 

(Sheet  1  of  2) 


92 


Load  Factor.  Nz  (g)  Pitch  Angle,  6  (radians) 


T-38A  USAF  $/N  63-81 35  QR03S  WEIGHT  -  1 0,378  LBS 

CRUISE  CONFIGURATION  CG:  1 6.5  %  MAC 
NO  EXTERNAL  STORES  YAW  AUGMENTATION:  OFF 


Time  (seconds) 


Time  (seconds) 


Figure  B3.  Time  History  Match,  Test  Point  2(a) 

(Sheet  2  of  2) 
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Pitch  Rate,  q  (radians  /  second)  Angle  of  Attack,  a  (radians) 


Figure  B4.  Time  History  Match,  Test  Point  2(b) 

(Sheet  1  of  2) 
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T-38AUSAFS/N  63-8135  GROSS  WEIGHT- 10,118  LBS 

CRUISE  CONFIGURATION  CG:  16.7  %  MAC 

NO  EXTERNAL  STORES  YAW  AUGMENTATION:  OFF 


Figure  B4.  Tine  History  Match,  Test  Point  2(b) 

(Sheet  2  of  2) 
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Pitch  Rate,  q  (radians  /  second)  Angle  of  Attack,  a  (radians) 


T-38AUSAF^63^135 
CRUISE  CONFIGURATION 
NO  EXTERNAL  STORES 


GROSS  WEIGHT  -  9,928  LBS 

CQ:  16.8  %  MAC 

YAW  AUGMENTATION:  OFF 


Figure  B5.  Time  History  Match,  Test  Point  2(c) 

(Sheet  1  of  2) 
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Pitch  Rate,  q  (radians  /  second)  Angle  of  Attack  a  (radians) 


T-38AUSAFS/N  63-8135 
CRUISE  CONFIGURATION 
NO  EXTERNAL  STORES 


GROSS  WEIGHT  -  11 ,41 8  LBS 

CG;18.1  %  MAC 

YAW  AUGMENTATION:  OFF 


Figure  B6.  Time  History  Match,  Test  Point  3(a) 

(Sheet  1  of  2) 


Load  Factor,  Nz  (g)  Pitch  Angle,  6  (radians) 


Figure  B6.  Time  History  Match,  Test  Point  3(a) 

(Sheet  2  of  2) 
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Pitch  Rate,  q  (radians  /  second)  Angie  of  Attack,  a  (radians) 


T-38A  USAF  S/N  63-81 35  GROSS  WEIGHT  -  1 0,91 8  LBS 

CRUISE  CONFIGURATION  CG:  1 7.3  %  MAC 

NO  EXTERNAL  STORES  YAW  AUGMENTATION:  OFF 


Time  (seconds) 


Figure  B7.  Time  History  Match,  Test  Point  3(b) 

(Sheet  1  of  2) 
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Load  Factor,  Nz  (g)  Pitch  Angle,  6  (radians) 


T-38A  USAF  S/N  63-8135 
CRUISE  CONFIGURATION 
NO  EXTERNAL  STORES 


GROSS  WEIGHT  -  10,918  LBS 

CG:  17.3  %  MAC 

YAW  AUGMENTATION:  OFF 


Figure  B7.  Time  History  Match,  Test  Point  3(b) 

(Sheet  2  of  2) 
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Pitch  Rate,  q  (radians  /  second)  Angle  of  Attack,  a  (radians) 


T-38AUSAFS/N  63-8135 
CRUISE  CONFIGURATION 
NO  EXTERNAL  STORES 


GROSS  WEIGHT  -  10,818  LBS 

CG:  17.1  %  MAC 

YAW  AUGMENTATION:  OFF 


Figure  B8.  Time  History  Hatch,  Test  Point  3(c) 

(Sheet  1  of  2) 
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T-38A  USAF  S/N  63-8135 
CRUISE  CONFIGURATION 
NO  EXTERNAL  STORES 


GROSS  WEIGHT  -  10,818  LBS 

CG:  17.1  %  MAC 

YAW  AUGMENTATION:  OFF 


Figure  B8.  Time  History  Natch,  Test  Point  3(c) 

(Sheet  2  of  2) 
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of  Attack,  a  (radians) 


Figure  B9.  Time  History  Match,  Test  Point  4(a) 

(Sheet  1  of  2) 
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Pitch  Angle.  6  (radians) 


Figure  B9.  Time  History  Match,  Test  Point  4(a) 

(Sheet  2  of  2) 
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Pitch  Rate,  q  (radians  /  second)  Angle  of  Attack,  a(radians) 


T-38A  USAF  S/N  63-6135 
CRUISE  CONFIGURATION 
NO  EXTERNAL  STORES 


GROSS  WEIGHT  -  1 0,91 8  LBS 

CG:  17.3  %  MAC 

YAW  AUGMENTATION:  OFF 


Figure  BIO.  Time  History  Match,  Test  Point  4(b) 

(Sheet  1  of  2} 


Figure  BIO.  Time  History  Match,  Teat  Point  4(b) 

(Sheet  2  of  2) 


107 


Pitch  Rate,  q  (radians  /  second)  Angle  of  Attack,  a  (radians) 


Figure  Bll.  Time  History  Match,  Test  Point  4(c) 

(Sheet  1  of  2) 
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Load  Factor,  Nz  (g)  Pitch  Angle,  6  (radians) 


T-38AUSAFS/N  63-8135 
CRUISE  CONFIGURATION 
NO  EXTERNAL  STORES 


GROSS  WEIGHT  -  10,618  LBS 

CG:  16.8  %  MAC 

YAW  AUGMENTATION:  OFF 


Tima  (seconds) 


Time  (seconds) 


Figure  Bll.  Time  History  Match,  Test  Point  4(c) 

(Sheet  2  of  2) 
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Phase  Angle  (degrees)  Bode  Magnitude  (dB) 


T-38A  USAF  S/N  63-81 35  AVERAGE  WEIGHT:  1 0.1 41  LBS 

CRUISE  CONFIGURATION  AVERAGE  CG:  1 6.7  %  MAC 
NO  EXTERNAL  STORES  YAW  AUGMENTATION:  OFF 


LEGEND 


Frequency  (radians  /  second) 


Figure  B12.  FRA  Results,  Test  Point  2(a) 
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Phase  Angle  (degrees)  Bode  Magnitude  (dB) 


T-38A  USAF  S/N  63-81 35  AVERAQ  E  WEIGHT:  1 1 ,051  LBS 

CRUISE  CONFIGURATION  AVERAGE  CG:  17.5  %  MAC 
NO  EXTERNAL  STORES  YAW  AUGMENTATION:  OFF 


LEGEND 


Frequency  (radians  /  second) 


Frequency  (radians  /  second) 


Figure  B13.  FRA  Results,  Test  Point  3(a) 
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Figure  B14.  FRA  Results,  Test  Point  4(a) 
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